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EXECUTIVE  SUMMARY 


There  are  more  than  cells  in  the  human  body,  and  every  second,  more  than  10  million  die 
and  are  replaced.  The  unity  of  an  organism  is  coordinated  by  many  levels  of  material,  chemical, 
and  electric  communication  between  the  cells.  Until  recently  progress  in  cell  biology  was 
limited  by  human  comprehension  of  complex  cellular  interactions  observed  from  laboratory 
experiments.  The  large  number  of  variables,  complex  functionality,  and  nonlinear  interactions 
in  cell  biology  events  severely  limit  human’s  capability  to  directly  comprehend  cellular 
behavior.  The  only  logical  way  to  compactly  describe  and  understand  such  complex  systems  is 
to  use  mathematics  and  computational  modeling.  Rapid  accumulation  of  genomic  and 
proteomic  data  from  novel  high  throughput  experimental  screening  technologies  demand  novel 
mathematical  approaches  and  models  to  process  and  interpret  massive  amounts  of  data.  These 
models  will  ultimately  address  a  variety  of  dynamic  intracellular  processes  ranging  from 
interactions  within  a  gene  regulation  network  to  intracellular  and  intercellular  signal 
transduction.  At  present,  however,  compact  elegant  mathematical  models  of  cellular  biology  are 
rare. 

In  the  last  two  decades,  computational  technologies,  such  as  Computational  Fluid  Dynamics 
(CFD)  and  Finite  Element  Method  (FEM)  based  Computational  Structures  Dynamics  (CSD), 
have  revolutionized  classical  engineering  disciplines  including  fluid  dynamics  and  structures 
mechanics.  We  believe  that  emerging  computational  technologies  are  poised  to  produce  a 
greater  revolution  in  medicine  and  biology.  Based  on  that  belief,  CFD  Research  Corporation 
(CFDRC)  formed  a  new  Computational  Medicine  and  Biology  (CMB)  division  in  2002  that  is 
determined  to  spearhead  that  revolution. 

Scientific  potentials  and  military  relevance  of  computational  biology  and  bioinformatics  have 
also  attracted  DARPA’s  interest.  The  DARPA  IPTO  ambitious  and  visionary  Bio-Computation 
(BioCOMP)  Program  to  develop  open  source  computational  framework,  called  Bio-SPICE, 
along  with  broad  range  of  cell  biology  modeling  tools  was  started  in  2001.  Bio-SPICE  stands  for 
the  Simulation  Program  for  Intra-Cellular  Processes.  The  BioCOMP  Program  was  tasked  to 
develop  novel  modeling  methods  and  simulation  software  for  broad  but  complementary  aspects 
of  cell  biology  such  as:  metabolic  pathways,  cellular  signaling,  cell  cycle  and  division, 
chemosensing,  chemotaxis,  cellular  colonies,  and  microorganisms.  The  CFDRC  role  in  the 
BioCOMP  Program  was  to  explore  multidimensional,  spatiotemporal  mathematical  modeling 
techniques  in  cellular  biology. 

Prior  to  this  project  CFDRC  was  involved  in  the  DARPA  BioCOMP  subcontract  to  the 
Lawrence  Berkeley  Laboratory  (LBL)  working  on  mathematical  formulation  and  computational 
exploration  of  spatiotemporal  cell  biology  models.  Summary  of  that  project  is  presented  in  the 
Appendix  A  of  this  report.  That  effort  provided  a  good  starting  point  for  the  present  project, 
whose  objective  was  to  develop  spatiotemporal  computational  cell/organ  biology  tools. 

The  focus  of  the  project  described  in  this  report,  was  to  develop  a  self-contained  software 
framework  solving  partial  differential  equations  governing  flow,  diffusion,  and  biochemistry 
applicable  for  spatiotemporal  modeling  of  cell  and  organ  biology  problems.  During  a  two-year 
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period  CFDRC  developed  Computational  Biology  (CoBi,)  software  tools  to  simulate  complex 
cell  and  organ  biology  problems.  The  code  was  written  in  modem  C++  programming  language 
using  state  of  the  art  numerical  techniques  for  solving  large  systems  of  partial  differential 
equations  (PDEs)  and  CFDRC  multiphysics  modeling  expertise.  A  graphical  user  interface 
(GUI)  for  CoBi,  JCoBi,  was  written  in  Java  and  interactive  3D  graphics.  CoBi  has  been 
designed  to  interact  with  the  other  Bio-SPICE  software  tools  using  Systems  Biology  Markup 
Language  (SBML,  http://sbml.org/index.psp  )  as  a  standard  interface. 

This  report  describes  the  mathematical  formulation  of  cellular  biophysics  and  biochemistry 
models  used  in  the  code  development,  numerical  methods  for  solving  large-scale  systems  of  PDE 
equations,  and  CoBi  software  stmcture.  The  code  has  been  successfully  applied  to  a  number  of 
cell  biology  problems.  This  report  describes  the  models  in  detail  and  presents  simulation  results 
for  several  cell  biology  problems  including:  bacterial  chemosensing  and  chemotaxis,  bacterial 
sporalation,  epidermal  growth  factor  receptor  (EGFR)  signal  transduction,  cellular  and  tissue 
calcium  oscillations,  cellular  and  tissue  oxygen  and  energy  metabolism,  morphogenesis  of  the 
yeast  cell,  and  perfusion  of  a  cell  in  an  organ.  The  CoBi  code  has  been  developed,  validated,  and 
demonstrated  on  several  cell  biology  and  organ  physiology  problems.  The  code  has  been  used  as 
the  starting  point  for  several  DARPA  and  DoD  projects  related  to  military  medicine,  DNA  based 
bio-computing,  personnel  protection,  and  biodefense. 

This  project  has  established  excellent  foundations  for  future  computational  medicine  and  biology 
projects  (CMB)  at  CFDRC  for  military  and  civilian  medical  applications.  After  two  years  of  this 
project,  CoBi  has  matured  to  the  point  that  it  can  be  used  for  solving  problems  not  only  in  cell 
biology  but  also  organ  biology  and  whole  body  (organism)  physiology.  CFDRC  has  invested 
significant  IR&D  resources  in  CoBi  code  development,  integration,  and  testing.  The  code  is 
being  used  as  a  starting  point  for  several  projects  with  DoD  focused  on  biotechnology, 
nanotechnology,  and  military  medicine. 
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1.0  INTRODUCTION 


1.1  Background 

CFD  Research  Corporation  (CFDRC)  has  been  involved  in  the  DARPA/IPTO  BAA  01-26 
BioCOMP  Program  between  January  2002  and  April  2006.  The  objective  of  the  DARPA 
BioCOMP  Program,  conceptualized  and  organized  by  Dr.  Sri  Kumar  of  DARPA/IPTO  was  to 
develop  a  computational  framework  that  enables  the  construction  of  sophisticated  models  of 
intracellular  processes  that  can  be  used  to  predict  and  control  the  behavior  of  living  cells.  One  of 
the  program  tracks  was  to  develop  Bio-SPICE  modeling  framework,  which  stands  for  the 
Simulation  Program  for  Intra-Cellular  Processes.  DARPA  has  initiated  the  project  to  develop 
computational  models  and  open  source  software  for  modeling  cell  biology.  Among  several 
projects,  typically  lead  by  US  academic  teams,  CFDRC’s  project  task  was  first  to  formulate  and 
explore  a  spatiotemporal  (multidimensional)  modeling  approaches  and  then  to  implement  and 
demonstrate  it  on  cell  and  organ  biology  problems. 

The  project  described  in  this  report  was  conducted  as  a  direct  contract  funded  by  DARPA  and 
coordinated  and  technically  directed  by  the  Air  Force  Research  Laboratory  (AFRL/IFTC  Rome 
NY).  Mr.  Clare  Thiem  of  AFRL/IFTC  was  the  DoD  Technical  Project  Manager  for  the  CFDRC 
BioCOMP  Program.  This  project  was  conducted  during  June  2004  to  May  2006.  The  overall 
objective  of  the  CFDRC  BioCOMP  Program  was  to  develop  a  self-contained  software 
framework  solving  partial  differential  equations  governing  flow,  diffusion,  and  biochemistry 
applicable  for  spatiotemporal  modeling  of  cell  and  organ  biology  problems.  The  specific 
objective  of  this  project  was  to  develop  Computational  Biology  (CoBi)  software  tools  and  a 
graphical  user  interface  (GUI)  using  state-of  the  art  numerical  techniques  and  novel 
programming  tools  (C++  and  Java)  to  simulate  complex  cell  and  organ  biology  problems.  CoBi 
has  been  designed  to  interact  with  the  other  Bio-SPICE  software  tools  using  Systems  Biology 
Markup  Language  (SBML,  httr)://sbml.org/index.r)sp  )  as  a  standard  interface.  The  prototype 
CoBi  software  developed  in  this  project  was  delivered  in  source  to  DARPA  (SRI  Bio-SPICE 
site)  in  February  2005  and  the  final  code  along  with  the  code  documentation  has  been  delivered 
to  AFRL/IFTC  at  the  end  of  the  project. 

1.2  Project  Objectives 

The  overall  objective  of  the  CFDRC  BioCOMP  Program  was  to  develop  a  software  toolkit  and 
modeling  framework  for  multidimensional,  spatiotemporal  comprehensive  modeling  of  cell 
biology. 

The  objective  of  the  present  follow  up  project  was  to  develop  a  self-contained  software 
framework  and  computational  cell  and  organ  biology  for  applications  in  military  medicine,  bio¬ 
computation,  biotechnology,  and  in  civilian  medical  and  biological  sciences.  The  experience 
gained  in  CFDRC’s  earlier  involvement  in  the  BioCOMP  Program,  described  in  Appendix  A,  in 
modeling  3D  cell  biology  problems  and  the  established  biomembrane  model  was  a  starting  point 
for  this  project.  CFDRC  proposed  to  adapt  the  existing  3D  model  of  cell  biomembrane  for 
complete  simulation  of  cell  biology  including  intracellular,  extracellular  and  membrane  bound 
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biophysical  and  biochemical  events.  The  CoBi  software  tools  should  be  able  to  solve 
generalized  adveetion-diffusion-reaetion  equations  governing  eell  biology  and  biophysies  using 
state-of-the-art  numerical  methods.  Thus,  the  objectives  for  this  project  were  as  follows: 

1 .  Develop  a  3D  unstructured  grid  solver  for  volumetric  diffusion-reaction  processes, 

2.  Integrate  the  existing  membrane  model  (developed  in  the  previous  DARPA  BioCOMP 
Program )  with  the  volumetrie  model  using  the  embedded  boundary  approach, 

3.  Develop  a  generalized  bioehemistry  solver  module  for  modeling  complete  cell  biology 
problems  including  extracellular,  intraeellular,  and  membrane  bound  bioehemieal  reaetions, 

4.  Establish  general  computational  biology  software  tools  for  modeling  eell,  colony,  tissue, 
and  organ  biology.  Provide  a  framework  for  3D  transient  solution  of  any  number  of 
adveetion-diffusion-reaetion  partial  differential  equations  in  eomplex  geometries  of  eellular 
and  eulture/tissue  structures, 

5.  Develop  interfaces  between  CoBi  and  other  Bio-SPICE  tools  and  eommunity  standard 
databases,  modeling  languages  (e.g.  SBML),  and  open  source  graphieal  post  processing 
tools  (e.g.  VizIT), 

6.  Develop  flexible  Graphieal  User  Interfaee  (GUI)  for  CoBi,  JCoBi,  to  enable  easy  model 
setup  and  simulation  eontrol, 

7.  Validate  the  software  on  several  benchmark  problems, 

8.  Demonstrate  CoBi  on  several  defense  related  eell  and  organ  biology  problems,  and  utilize  it 
for  other  DoD  projeets, 

9.  Document  CoBi  and  deliver  it  to  DARPA  and  to  DoD  Laboratories. 

It  is  felt  that  all  of  the  major  objeetives  of  this  projeet  have  been  suecessfully  accomplished.  The 
CoBi  code  has  been  developed,  validated,  and  demonstrated  on  several  cell  biology  and  organ 
physiology  problems.  The  code  has  been  used  as  the  starting  point  for  several  DARPA  and  DoD 
projeets.  CFDRC  is  using  CoBi  as  the  enabling  teehnology  for  projeets  related  to  military 
medieine,  DNA  based  biocomputing,  personnel  proteetion,  and  biodefense.  The  eode  has  been 
delivered  to  interested  DoD  Laboratories.  The  following  sections  describe  the  model,  software, 
and  examples  from  comprehensive  simulation  studies. 
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2.0 


COMPUTATIONAL  CELLULAR  BIOLOGY 


2.1  Relevant  Introduction  to  Cellular  Biology 

Living  systems,  from  simple  organisms  to  the  human  body,  eontain  organs.  Organs  are 
composed  of  tissues,  tissues  consist  of  cells,  and  cells  are  formed  of  organelles,  and  molecules. 
Figure  2.1-  the  biological  cell  is  the  basic  unit  of  all  living  organisms  (for  detailed  discussions 
of  cell  biology  see  e.g.  Lodish  2000).  There  are  more  than  300  trillion  cells  in  the  human  body, 
and  every  second  of  every  day,  more  than  10  million  die  and  are  replaced.  The  unity  of  an 
organism  is  coordinated  by  many  levels  of  material,  chemical,  and  electric  compunctions. 
Specific  organ  tissues  and  individual  cells  are  surrounded  by  molecular  membranes  produced  by 
cells.  The  biological  universe  consists  of  two  types  of  cells,  shown  schematically  in  Figure  2.1 
prokaryotic  cells,  and  eukaryotic  cells. 

Prokaryotes,  such  Escherichia  Coli  bacteria,  are  single  cell  organisms,  which  have  simple 
internal  organization  lack  a  nucleus,  and  the  genetic  material  (DNA)  is  directly  suspended  inside 
the  cell  cytoplasm.  Eukaryotes,  comprising  plant  and  animal  kingdoms,  have  complex  internal 
structure  with  various  organelles  and  membrane  limited  nucleus  that  is  kept  separate  from  the 
cytoplasm  by  a  double  membrane  structure.  The  cytoplasm  contains  the  rest  of  the  organelles 
and  the  area  of  the  cytoplasm  outside  of  the  individual  organelles  is  called  the  cytosol.  Each 
organelle,  surrounded  by  a  membrane,  contains  a  collection  of  specific  enzymes  and  plays  a 
unique  role  in  the  growth  and  metabolism  of  the  cell.  The  most  important  organelles  are: 
nucleus  -  containing  DNA  and  other  replication  molecules,  mitochondria  -  in  which  most  of  he 
cell  energy  metabolism  takes  place,  endoplasmic  reticula  -  build  of  complex  membranes  within 
which  glycoproteins  and  lipids  are  synthesized,  Golgi  vesicles  -  which  direct  membrane 
constituents  to  appropriate  places  in  cells,  peroxisomes  -  in  which  fatty  acids  and  aminoacids 
are  degraded,  and  lysosomes  -  which  also  degrade  internal  and  foreign  molecules. 


Figure  2.1:  Examples  of  mammalian  (eukaryotic)  and  bacterial  (prokaryotic)  cells  (with  NIH 
permission),  htW://www. ncbi. nlm.nih. sov/About/primer/senetics  cell. html. 

Biological  cells  are  surrounded  by  a  membrane,  phospholipid  bilayer,  designed  to  keep  the  cell 
interior  together  and  protect  the  cell  from  the  surrounding  environment.  The  biochemical 
processes  of  cellular  life  takes  place  in  the  cytosol  and  cells  communicate  with  the  surrounding 
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and  neighbor  cells  via  reactions  on  the  membrane  and  via  the  trans-  membrane  molecule 
exchange.  The  cytosol  is  the  largest  structure  in  the  cell.  It  composes  54%  of  the  cells  total 
volume.  The  cytosol  contains  thousands  of  enzymes  that  are  responsible  for  the  catalyzation  of 
glycolysis  and  gluconeogenesis  and  for  the  biosynthesis  of  sugars,  fatty  acids,  and  amino  acids. 
The  cytosol  takes  molecules  and  breaks  them  down,  so  that  the  individual  organelles  can  use 
them. 

The  cytosol  contains  a  skeletal  structure,  called  the  cytoskeleton.  Cytoskeleton  gives  the  cell  its 
shape  and  allows  it  to  organize  many  of  the  chemical  reactions  that  occur  in  the  cytoplasm. 
Additionally,  the  cytoskeleton  can  aid  in  the  movement  (motility)  of  the  cell.  It  consists  of 
protein  fdaments:  actin  filaments  and  the  microtubules.  The  actin  is  responsible  for  contraction 
(like  in  muscles)  and  the  microtubules  are  for  structural  strength  of  the  cell.  Cytoskeletal  fibers 
also  control  movement  of  structures  within  the  cell. 

To  survive  cells  have  to  carry  out  numerous  functions  ranging  from  replication  and  energy 
conversion  to  molecule  transport  and  the  various  complicated  cascades  of  biochemical  reactions, 
signaling,  used  in  cellular  communication.  Some  of  the  most  important  functions  are:  Molecular 
transport,  DNA  replication.  Reproduction,  Protein  synthesis,  and  Metabolism  and  Signaling. 

One  of  the  primary  goals  of  all  living  organisms  is  to  survive  and  reproduce  -  typically  by  cell 
division.  Cell  division  occurs  rapidly  in  living  organisms.  For  example,  in  an  adult  human, 
millions  of  cells  divide  each  second  to  maintain  homeostasis  (the  proper  balance  in  cells).  Cells 
can  reproduce  in  two  ways,  mitosis,  and  meiosis.  In  mitosis,  the  resulting  daughter  cell  is  an 
identical  clone  of  the  original  cell.  Mitosis  is  mostly  used  by  somatic  cells  (cells  of  the  body). 
Meiosis,  however,  is  a  form  of  sexual  reproduction  and  only  occurs  in  reproductive  cells. 

Mitotic  cell  division  is  an  ordered  set  of  events,  cell  cycle,  culminating  in  cell  growth  and 
division  into  two  daughter  cells.  It  begins  with  interphase,  when  the  cell  replicates  all  of  its 
genomic  and  cytoplasmic  material  and  prepares  for  division.  After  preparation  is  complete,  the 
cell  enters  the  4-phased  mitosis.  In  mitosis,  the  cell  sequentially  goes  through  prophase, 
metaphase,  anaphase,  and  telophase.  Immediately  after  the  completion  of  telophase, 
cytokenesis  is  initiated  to  end  the  cell  division  by  literally  separating  the  cell  in  two.  Figure 
2.2(a)  schematically  illustrates  the  cell  division  and  cell  cycle  processes. 

The  cell  cycle  stages.  Fig  2.2(b),  are  G1-S-G2-M.  The  G1  stage  stands  for  "GAP  1".  The  S 
stage  stands  for  "Synthesis".  This  is  the  stage  in  which  DNA  replication  occurs.  The  G2  stage 
stands  for  "GAP  2".  The  M  stage  stands  for  "mitosis",  and  is  when  nuclear  (chromosomes 
separate)  and  cytoplasmic  (cytokinesis)  divisions  occur. 


Figure  2.2(a)  Figure  2.2(b) 

Figure  2.2:  Schematics  of  events  during  cell  division  and  cell  cycle  (with  NIH  permission, 
httv://www. ncbi. nlm. nih. 20v/About/vrimer/2enetics  cell. html . 
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Cytokinesis  is  the  final  event  of  the  cell  cycle  and  is  the  process  that  divides  one  cell  into  two 
daughter  cells.  Figure  2.3  illustrates  how  different  organisms  conduct  cytokinesis.  In  animal 
cells,  the  division  site  is  first  chosen,  generally  at  the  cell  equator,  and  subsequently  the  cleavage 
furrow  is  assembled  at  the  division  site  (Figure  2. 3D).  The  furrow  contains  actin,  myosin,  and 
other  proteins  that  are  organized  into  a  contractile  ring  called  the  actomyosin  ring.  The  ring  then 
ingresses  or  contracts,  generating  a  membrane  barrier  between  the  cytoplasmic  contents  of  each 
daughter  cell.  The  ingressing  furrow  constricts  components  of  the  spindle  midzone  into  a 
focused  structure  called  the  midbody.  In  the  final  cytokinetic  event,  called  abscission,  the  furrow 
“seals”,  generating  two  completely  separate  cells. 


B 
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Figure  2.3:  General  mechanisms  of  cytokinesis  in  eukaryotes.  (B)  In  buddins  yeast  cells,  the 
ring  is  positioned  at  the  interface  between  the  mother  cell  and  daughter  bud,  (C)  In  fission  yeast 
and  animal  cells,  the  contractile  ring  is  centrally  placed,  as  both  cell  types  divide  by  medial 
fission.  Both  budding  and  fission  yeasts  synthesize  a  division  septum,  which  eventually 
degrades,  resulting  in  physical  cell  separation.  (D)  In  animals,  the  ingressing  furrow  constricts 
the  spindle  midzone  components  into  a  dense  structure-  midbody.  Guertin2002.  (with  NIH 
permission),  http://pubmedcentral.com/articlerender.fcsi?artid= 120788 . 


Live  biological  cells  exhibit  an  amazing  capability  to  sense  the  surrounding  environment  and 
ability  to  respond  to  it.  By  rearranging  their  interior  (cytoskeleton,  organelles,  membrane) 
mammalian  cells  grow,  divide,  move,  expand/contract,  and  much  more.  Cell  movement  toward 
(or  away)  from  an  external  stimulus,  cell  taxis,  is  induced  by  cell  sensory  organs,  which  activate 
cell  asymmetry,  deformation  and  motility  mechanisms.  Cell  taxis  stimulated  by  chemical  agents 
-  chemotaxis  allows  them  to  detect  the  direction  (gradient)  and  proximity  of  an  external 
chemical  stimulus  such  as  nutrient  and  energy,  [Bray  2001]. 


Within  a  complex  organism  organs  and  individual  cells  communicate  by  exchanging  chemical 
messengers  traveling  from  one  cell  to  another  by  direct  contact  or  through  the  extracellular 
space  (interstitium).  The  chemical  signals  received  by  target  cells  respond  by  activating 
signaling  pathways  involving  intracellular  proteins  (metabolic  enzymes,  gene  regulatory 
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proteins,  ion  channels,  or  cytoskeletal  proteins).  Signaling  molecules,  also  called  ligands 
(protein,  small  peptide,  amino  acid,  nucleotide,  steroid,  hormone)  are  secreted  from  one  cell  in 
response  to  a  specific  stimulus  travel  to  a  target  cell,  where  they  bind  to  specific  receptors  on 
their  membrane  and  elicit  an  intracellular  response. 

Within  the  cell  the  main  method  of  signal  transduction  occurs  through  structural  changes  of 
pathway  components.  A  given  protein  will  affect  the  conformation  of  one  or  several  other 
proteins,  activating  or  inhibiting  those  proteins  and  thus  propagating  the  signal  down  the 
pathway.  The  trigger  for  signal  propagation  often  occurs  with  the  binding  of  the  signaling 
molecule  to  the  receptor,  which  causes  a  conformational  change  in  the  receptor.  Subsequently, 
cytosolic  regions  of  the  receptor  are  activated,  making  them  active  targets  for  intracellular  and 
membrane  associated  proteins.  The  biochemical  signal  transduction  in  the  cell  is  a  complex  and 
combinatorial  process  involving  switches,  integrating  centers,  feedback  loops  and  crosstalk 
between  pathways.  Figure  2.4  presents  example  schematics  of  cell  signaling  pathways  analyzed 
at  CFDRC  for  modeling  cellular  neurotoxicity  problems.  Obtaining  a  detailed  understanding  of 
the  dynamics  of  a  biochemical  reaction  is  a  formidable  challenge. 


Figure  2.4:  Schematic  of  organophosphate  neuron  toxicity  cellular  signaling  pathways,  from 
CDFRC  dose-dependent,  analysis  of  signaling  and  activation,  neurotransmitter  release,  and 
metabolic  response  in  neurotypical  PCI 2  cells,  Jenkins  J.  and  Hood  J.  2006. 
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2.2  Biological  and  Computational  Challenges  in  Cellular  Biology 

Rapid  progress  in  experimental  biology  and  medicine,  particularly  at  molecular  and  cellular 
levels,  is  generating  an  enormous  amount  of  “bio-informatics”  data.  The  raw  data  in  digital  or 
visual  imaging  form  is  analyzed  and  interpreted  by  highly  skilled  experts  who  formulate 
“models”  of  cellular  events.  These  “models”  are  typically  expressed  in  a  long  narrative  form 
and  described  in  biomedical  journals.  In  molecular  and  cellular  biology  compact  elegant 
mathematics  based  models  are  rare. 


In  the  last  few  years  it  has  been  recognized  that  further  understanding  of  biology  will  be  possible 
only  with  a  more  comprehensive  multidisciplinary  approach  involving:  mathematics, 
informatics,  computational  biophysics,  in  addition  to  more  traditional  disciplines  such  as 
biochemistry,  biology,  physiology,  and  pathology,  etc.  Joel  E.  Cohen  of  the  Rockefeller  & 
Columbia  Universities  elegantly  projected  progress  in  biology  as: 

"^Mathematics  Is  Biology ’s  Next  Microscope,  Only  Better; 

Biology  Is  Mathematics  ’  Next  Physics,  Only  Better  ” 

A  Computational  and  bioinformatics  approach  to  biology  and  medicine  has  enormous  potential 
that  has  not  yet  been  completely  revealed.  In  the  last  two  decades,  computational  technology 
(hardware  and  software)  has  been  revolutionized  and  classical  engineering  disciplines  such  as 
fluid  dynamics  or  structures  mechanics  evolved  into  Computational  Fluid  Dynamics  (CFD)  and 
Finite  Element  Method  (FEM)  based  Computational  Structures  Dynamics  (CSD).  New 
industries  have  emerged  and  CFDRC  is  one  of  the  best  examples  of  US  based  commercial 
ingenuity  and  business  successes.  It  is  believed  that  similar  scientific  and  business  revolutions 
are  approaching  in  Computational  Medicine  and  Biology  (CMB).  In  2002  CFDRC  formed  a 
new  CMB  team  determined  to  spearhead  that  revolution. 

Scientific  potentials  and  military  relevance  of  computational  biology  and  bioinformatics  have 
also  attracted  the  interest  of  DARPA.  The  DARPA  IPTO  ambitious  and  visionary  BioCOMP 
Program  to  develop  an  open  source  computational  framework  and  modeling  tools  for  cell 
biology  started  in  2001.  This  framework  is  called  Bio-SPICE  which  stands  for  the  Simulation 
Program  for  Intra-Cellular  Processes.  The  BioCOMP  Program  involves  experts  from  a  broad 
range  of  disciplines,  including  biology,  biochemistry,  biophysics,  computer  science, 
mathematics,  and  bioengineering.  Several  teams  were  tasked  to  develop  novel  modeling 
methods  and  simulation  software  for  broad  but  complementary  aspects  of  cell  biology  such  as: 
metabolic  pathways,  cellular  signaling,  cell  cycle  and  division,  chemosensing,  chemotaxis, 
cellular  colonies,  and  microorganisms.  New  standards  of  model  description  languages  (System 
Biology  Modeling  Language,  SBML)  were  being  pursued  to  enable  model  exchange,  reuse  and 
archiving. 

In  the  current  BioCOMP  Program  almost  all  Bio-SPICE  development  activities  concentrate  on 
the  development  of  biochemistry  (cell  cycle,  signaling,  metabolism,  genomics,  etc)  or 
bioinformatics  (e.g.  microarray  data  processing)  use  Ordinary  Differential  Equations  (ODEs) 
and  as  a  consequence  describe  cells  as  a  “well  stirred  reactor”.  Spatial  multidimensional  cell 
modeling  using  Partial  Differential  Equations  (PDEs)  and  accounting  for  cell  morphology  are 
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much  more  challenging  and  not  well  represented  in  the  BioCOMP  Program.  Representing  cells 
as  “well  stirred  reactors”,  used  in  all  Bio-SPICE  biochemistry  modeling  tools,  is  clearly 
inadequate  even  for  the  simplest  cells.  The  reaction  mechanisms  and  rate  constants  in  those 
models  are  “tuned”  based  on  empirical  data  to  fit  experiments.  Those  rate  constants  are 
“contaminated”  by  other  biophysics  processes  not  accounted  for  explicitly  in  the  models  such  as: 
diffusion,  advection,  directed  transport  (e.g.  by  molecular  motors),  cytoplasm  cytoskeleton  and 
membrane  mechanics,  electrostatics,  and  others.  Such  models  are  also  not  scalable,  rate 
constants  may  be  time  dependent,  or  additional  stochastic  processes  are  added  to  cover  for 
missed  biophysics.  It  is  likely  that  all  the  biochemistry  reaction  mechanisms  rate  constants  will 
have  to  be  recalibrated  once  used  in  more  complete  spatiotemporal  cellular  biology  models. 

Therefore  there  is  an  urgent  need  for  prompt  development  of  spatiotemporal  cell  biology 
models.  The  objective  of  the  CFDRC  project  was  to  develop  a  software  toolkit  for 
multidimensional,  spatiotemporal  simulation  of  cell  biology. 

The  BioCOMP  Program  will  have  enormous  implications  not  only  in  basic  scientific  research 
but  also  more  importantly  in  the  US  national  commercial  business  and  defense.  Validated 
computational  cell  biology  modeling  tools  will  be  used  to  assess  novel  drugs,  novel  therapies, 
safety,  toxicity,  and  efficacy  of  drugs,  drugs  designed  for  individual  patients,  etc.  From  the 
military  perspective,  better  biochemical  sensors,  better  vaccines  for  bio-agents,  better  treatment 
of  trauma  of  injured  soldiers,  etc.  can  be  designed.  Such  tools  may  also  influence  paradigms  of 
biocomputing. 

2.3  Assessment  of  Cell  Biology  Spatial  Modeling  Methods 

Most  outstanding  progress  has  been  achieved  in  computational  cellular  biochemistry  where  cell 
biological  pathway  models  are  described  by  large  systems  of  ODEs  [Fall  2000].  Public  domain 
software  tools  are  being  developed  e.g.  by  Arkin  (1998)  (Bio-SPICE  http://BioSPICE.lbl.gov/)  in 
US  and  by  Tomita  (E-Cell,  http://www.e-cell.org)  see  e.g.  Tomita  [2001],  and  Takahashi  [2003] 
in  Japan.  Levchenko  2002  has  recently  presented  a  comprehensive  example  of  an  ODE  based 
model  of  cellular  chemotaxis.  However  they  all  lack  spatial  dimension  and  are  at  best 
“compartmentalized”. 

Multidimensional  models  of  cells  have  been  developed  by  Loew  [Schaff  1999,  Fink,  2000, 
Slepchenko  2003]  and  termed  a  Virtual  Cell.  ID  and  2D  spatiotemporal  computational  models 
of  chemotaxis  have  also  been  reported  by:  Bottino  [2000,  2002],  the  Lauffenburger  group  at  MIT 
[Narang  2001],  Painter  [2000],  Mogilner  [2002],  and  Othmer  [2002].  The  Virtual  Cell  model 
provides  a  formal  framework  for  modeling  biochemical  and  biophysical  and  transport 
phenomena  in  2D  and  probably  3D  (never  published)  cell  geometries  while  considering  the 
subcellular  localization  of  the  major  cellular  organelles  (nucleus,  endoplasmic  reticulum).  Cell 
geometry  can  be  imported  for  microscopy  images.  Geometrical  shapes  of  the  cell  have  been 
represented  with  the  Cartesian  grid  and  curved  shapes  of  the  membrane  have  to  be  represented  as 
“stair  steps”.  Figure  2.5  presents  an  example  of  the  Virtual  Cell  simulation  results,  the  IP3- 
mediated  calcium  dynamics  in  differentiated  NlE-1 15  mouse  neuroblastoma  cells  with  complex 
neuronal  morphologies  [Fink  2000].  The  results  point  to  an  important  role  for  cellular  geometry 
in  controlling  the  spatial  and  temporal  patterns  of  intracellular  signals. 
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Spatiotemporal  models  of  biological  cells  have  also  been  applied  to  cellular  chemosensing  and 
chemotaxis.  Simulation  results  have  been  demonstrated  for  the  2D  mechano-chemical  cell 
motion  model,  quasi  2D  chemical  chemoattractant  gradient  sensing  model,  and  detailed 
equations  for  nucleation,  growth/capping  by  polymerization,  and  diffusion  of  actin  filaments  in 
cells.  For  physiological  parameters  the  models  show  remarkable  agreement  with  experimental 
observations.  A  2D  PDE  “fluid  based”  chemotaxis  model  has  been  recently  developed  by 
Bottino  2000,  2002,  where  reaction,  diffusion,  advection  equations  were  solved  on  triangular  and 
Voronoi  polyhedral  (2D)  meshes.  Figure  2.5  illustrates  the  CFDRC  simulation  of  in  vitro 
cellular  chemotaxis  events.  A  cellular  culture  dish  erythrocytes,  a  white  blood  cell,  and  a  motile 
bacterium  was  setup  to  simulate  the  relative  motion  of  the  bacterium  and  the  neutrophil  which  is 
stimulated  by  the  chemokines  released  by  the  bacterium,  A  Lagrangian  -Eulerian  model  has  been 
used  to  solve  movement  of  cells  and  diffusion  of  a  chemokine  in  the  culture. 


Figure  2.5:  CFDRC  simulation  of  a  bacterial  and  neutrophil  chemotaxis  in  an  in  vitro  cell 

culture. 
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3.0  SPATIOTEMPORAL  MODEL  OF  CELLULAR  BIOLOGY 


3.1  Computational  Cell  Biology 

Substantial  progresses  over  the  past  three  decades  in  biochemistry  and  molecular  cell  biology, 
coupled  with  emerging  high  throughput  techniques  for  detecting  protein-protein  interaction,  have 
opened  a  new  era  in  signal  transduction  research.  At  the  same  time  however,  recent  analyses  of 
cellular  signaling  pathways  have  revealed  that  they  are  extremely  complex  nonlinear  molecular 
processes  occurring  throughout  the  cell  volume  and  involve  multiple  species  and  signaling 
cascades.  Because  of  their  size  and  complexity,  these  networks  are  often  too  complicated  for  the 
human  mind  to  organize  and  analyze.  It  has  become  necessary  to  develop  mathematical  models 
to  understand  the  behavior  of  signaling  networks  controlling  cellular  biology. 

The  use  of  mathematical  models  and  computer  simulation  has  been  successfully  used  in 
pharmacology  and  drug  discovery  to  study  ligand-receptor  interactions, 
pharmacokinetics,  and  pharmacodynamics.  Computational  modeling  of  cellular  metabolic  and 
signaling  pathways,  the  main  focus  of  the  BioCOMP  Program,  has  been  explored  using  four 
major  modeling  frameworks:  chemical  kinetics  models,  compartmental  models,  stochastic 
models,  and  diffusion-reaction  models.  The  most  familiar,  simplest,  and  most  commonly 
followed  (in  the  BioCOMP  Projects)  modeling  framework  is  the  chemical  kinetics  approach.  In 
this  approach,  the  cell  is  treated  as  a  “well-stirred  reactor”.  The  dynamics  of  a  cellular  signaling 
in  that  approach  is  described  by  a  system  of  Ordinary  Differential  Equations  (ODEs)  in  the  form: 


- =  Generation  -  Consumption  (3-1) 

dt 


Where  C  is  the  species  concentration  and  the  right  hand  side  (RHS)  typically  involves  several 
terms  in  which  specie  C  participates  arranged  in  two  groups:  generation  of  specie  C  and  its 
consumption.  The  generation  and  consumption  term  can  be  a  constant  (e.g.,  synthesis),  first 
order  reactions  (e.g.,  degradation),  or  nonlinear  (e.g.,  second  order  reactions  or  Michelis-Menten 
kinetics  for  enzymes).  For  large  signaling  networks  such  as  those  shown  in  Figure  2.4  and 
Figure  2.5,  in  the  preceding  chapter,  a  large  system  of  ODEs,  each  having  large  number  of 
generation  and  consumption  term  on  the  RHS.  In  the  last  few  years,  advanced  numerical 
techniques  and  robust  software  tools  have  been  developed  for  solving  large  systems  of  ODEs 
such  as  those  resulting  from  biochemical  signaling  pathways.  Several  of  them  have  been 
developed  and  adapted  for  modeling  cell  biology  problems  during  the  DARPA  Bio-SPICE 
program  e.g.  JigCELL  software  developed  at  Virginia  Tech  by  Professor  John  Tyson  group. 

To  understand  the  cellar  biochemical  dynamics  that  arises  from  variation  in  both  space  and  time, 
either  compartmental  or  partial  differential  equation  (PDE)-based  models  are  needed.  Like 
chemical  kinetics  models,  compartmental  models  are  based  on  ODEs.  However,  instead  of 
treating  a  cell  as  a  single  well  stirred  reactor,  compartmental  models  divide  the  cell  into 
compartments  such  as:  membrane,  cytoplasm,  mitochondrion,  nucleus,  or  extracellular  space), 
collect  the  equation  describing  biochemical  kinetics  of  species  in  individual  compartments  and 
allow  species  exchange  between  compartments. 
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Compartmental  modeling  has  been  used  extensively  in  modeling  the  network  involved  in  protein 
secretion  and  trafficking,  and  developmental  biology.  For  the  cases  where  explicit  dependence 
on  the  spatial  variables  is  needed,  a  partial  differential  equation  (PDE)-based  model  is  needed. 
The  spatiotemporal  dynamics  of  cell  biochemical  signaling  can  be  described  by  the  general 
advection-diffusion-reaction  equation: 

dC 

^  +V-(UC)=  yDVC+  R  (3.2) 

^  '  Convection  '  '  Diffusion  '  Reaction  Source 

Transient 


where,  C  is  the  specie  concentration,  U  is  the  fluid  flow  velocity,  D  is  the  diffusivity  of  the 
specie  C  in  the  cellular  medium  and  R  is  the  generalized  source  term.  As  in  the  simple  kinetics, 
models  the  source  term  may  consist  of  several  generation  and  consumption  terms.  Note  that 
equation  (3.2)  includes  several  terms  for:  transient,  convection,  diffusion,  and  source/sink  term 
due  to  chemical  reactions.  Equation  (3.2)  is  the  PDE  because  the  concentration  depends  on  three 
spatial  directions  (x,  y,  z  coordinates)  and  on  time.  Unlike  in  the  ODE  equations,  PDEs  are 
solved  on  a  computational  grid.  Generation  of  3D  computational  grid  for  a  biological  cells  is  a 
challenging  problem  in  itself  as  the  cell  geometry  may  be  changing  in  time  and  may  be 
influenced  by  the  cell  biological  events  e.g.  chemotaxis.  In  such  cases,  cell  geometry  and  mesh 
need  to  be  computed  as  part  of  the  cell  biology  problem.  Numerical  solutions  of  even  large  scale 
ODE  systems  are  very  fast  (seconds  to  minutes  on  powerful  PCs).  Numerical  solutions  of  large 
scale  time  dependent  3D  PDEs  is  much  more  computationally  intensive.  Depending  on  the 
problem,  it  may  take  only  few  minutes  or  as  long  as  several  days  or  weeks.  There  are  very  few 
spatiotemporal  cell  biology  tools,  one  is  the  Virtual  Cell  (http://www.nrcam.uchc. edu/1  and  the 
other  CoBi  tool  developed  in  this  project. 

3.2  Basic  Transport  Equations 

Spatiotemporal  computational  modeling  of  flow,  diffusion,  and  reaction  processes  in  complex 
geometries  requires  several  preparatory  steps.  The  first  and  most  important  is  to  define  the 
objective  of  the  task  and  the  expected  result.  A  step-by-step  approach  is  strongly  recommended 
by  first  solving  simplified  but  related  problems  and  adding  complexity  after  each  successful  step. 
A  typical  modeling  procedure  involves  the  following  steps: 

1 .  Create  or  import  the  model  geometry  and  computational  grid, 

2.  Select  the  models  (e.g.  flow,  species  transport,  biochemistry  pathways,  etc), 

3.  Specify  volume  conditions,  boundary  conditions,  initial  conditions, 

4.  Specify  material  properties,  and  other  physical  constants, 

5.  Specify  the  solution  control  parameters  e.g.  number  of  times  steps,  under-relaxation, 

6.  Calculate  a  solution  using  CoBi. 

7.  Analyze  the  results  using  visualization  tools  e.g.  Visit  tools  from  LBL, 

After  step  7,  you  may  want  to  return  to  the  previous  steps  refining  mesh,  models,  or  parameters. 

All  models  in  CoBi  are  formulated  using  an  SI  unit  system  which  will  be  used  in  the  description 
of  the  mathematical  models.  However  using  the  appropriate  model  setup  the  user  can  change  the 
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system  and  even  use  a  nondimensional  formulation  -  this  requires  a  good  understanding  of  all 
terms  of  the  transport  equations,  models,  and  properties  used.  If  you  read  an  external  grid  file 
with  X,  Y,  Z  coordinates  of  grid  vertices  we  assume  it  is  setup  in  meters.  If  not,  one  will  need  to 
scale  the  grid  using  the  preprocessor. 

Throughout  this  section,  for  compactness,  we  will  use  the  vectorial  notation  to  describe  partial 
differential  equations.  The  reader  will  notice  that  most  of  them  have  a  very  similar  form  as 
shown  below: 


Transient 


+y-{u-c)  = 

V 

Convection 


V-Z)VC+  S 

^ ^  ^ 
Diffusion  Source 


(3.3) 


where  C  is  the  transported  property  e.g.,  a  concentration,  U  is  the  convective  velocity  vector,  D 
is  the  diffusion  coefficient,  and  S  is  the  source  term  e.g.  chemical  reaction  rate  representing 
formation  and  consumption  of  species  C.  The  nabla  operator,  V ,  for  a  scalar  variable  C  is  called 
gradient  for  a  vector  variable  U  is  called  divergence: 


VC  =  v-u  = 

Gradient  5xy  5z  Divergence  5xy  dz 


(3.4) 


CoBi  has  several  unique  modeling  components  related  to  cell  and  tissue  biology.  One  of  the 
most  important  is  the  membrane  model.  The  membrane  model  and  its  interaction  with  the 
volume  models  of  the  cytoplasm  and  extracellular  space  will  be  described  in  a  separate  section  in 
this  chapter. 

Spatiotemporal  computational  cell  biology  involves  a  coupled  solution  of  cell  biophysics  and 
biochemistry.  Cell  biophysics  may  involve  several  disciplines  including  fluid  flow,  heat  transfer, 
structural  mechanics,  electrostatics,  and  others.  They  are  all  mathematically  described  by 
complex  PDEs.  The  biochemistry  is  described  by  a  system  of  species  transport  equations. 
Equation  3.2. 

The  convective  diffusive  transport  in  cells  and  tissues  occurs  in  a  non-homogeneous  porous 
media.  The  governing  equations  of  fluid  (e.g.  water)  flow  in  the  biological  cells  and  tissues  are 
the  Navier-Stokes  equations,  modified  for  a  porous  medium,  and  consist  of  the  fluid  continuity 
and  momentum  equations.  Transport  of  multiple  chemical  species  such  as  signaling  molecules 
or  drugs  is  described  using  convective-diffusive  transport  equations  later  in  this  document.  The 
fluid  transport  equations  are  given  below: 

Continuity  Equation 

^(ep)  +  V*(€pV)  =  0  (3.5) 

ot 


In  this  equation,  the  first  term  is  the  mass  change  (density  or  volume)  in  the  control  volume.  It  is 
zero  for  incompressible  flow  and  if  the  computational  mesh,  is  not  moving  or  deforming.  The 
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flow  of  water  under  normal  pressure  and  temperature  is  a  good  example  of  incompressible  flow. 
The  second  term  is  the  mass  flux  across  surrounding  surface  of  the  same  control  volume. 

Fluid  Momentum  Equations 


— (£-/?U)  +  V  •  (i-pUU)  =  -sVp  +  V  •  (sx)  +  ^-f^  (3.6) 

dt  K 

The  first  term  on  the  left  is  the  rate  of  change  of  momentum;  the  second  term  is  the  convection 
term.  The  first  term  on  the  right  side  is  the  pressure  gradient  term,  the  second  term  is  viscous 
shear  stress  term,  and  the  last  term  is  the  resistance  term  arising  due  to  the  porous  medium.  In 
above  equations  8  is  the  porosity  (1.0  open  flow,  0.0  fully  blocked  volume),  p  is  the  density,  p  is 
the  pressure,  U  is  the  water  perfusion  velocity  vector,  x  is  the  shear  stress,  and  k  is  permeability. 
A  “continuum  approximation”  is  assumed  for  the  cellular  space  and  organ  tissue 
approximating  the  cells  and  intercellular  space  using  the  porosity  concept. 

Cellular  biology  involves  complex  interactions  between  the  biophysical  and  biochemical 
processes  occurring  in  the  extracellular  space,  on  cell  membranes,  and  within  intracellular  space. 
Vast  number  of  species  (metabolic  nutrients  and  product,  signaling  molecules,  drugs)  are 
created,  transported,  and  processes  in  various  regions.  Computational  modeling  of  these 
phenomena  involves  the  solution  of  species  transport/balance  equations  in  extracellular, 
membranes  and  intracellular  spaces  within  the  tissue.  The  general  transport  equation  for  species 
in  perfused  tissue  is  described  by  a  partial  differential  equation  (PDE),  Equation  (3.7),  and 
includes  several  terms  for:  transient,  convection,  diffusion,  and  source/sink  term  due  to  chemical 
reactions  and  can  be  written  as  follows: 


dsC 

Transient 


+  V  •  (^UC)  = 

Convection 


y  ■  sDS/C  + 

Diffusion  Reaction  Source 


(3.7) 


where,  U  is  the  perfusion  velocity,  D  is  the  diffusivity  of  drug  in  the  medium  and  s  is  the 
porosity  of  the  medium. 

The  species  diffusion  coefficient  D  in  equation  3.7  represents  “diffusivity  in  media”.  The 
experimentally  measured  total  tissue  effective  diffusivity  will  need  to  be  scaled  by  porosity  D/s. 
The  species  diffusivity  value  depends  on  several  parameters  representing  intercellular  space  and 
the  transported  species. 

Typically,  the  transport  and  biochemistry  terms  in  the  species  balance  PDE  equations  are  solved 
using  an  operator  splitting  technique  by: 

1.  First  solving  the  system  of  ODE  equations  composed  of  transient-reaction  terms  for  all 
species,  and  then 

2.  Solving  all  terms  of  the  PDE  equations  with  source  terms  computed  from  step  1 . 
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An  important  part  of  the  spatial  cell  modeling  is  the  geometry  definition  and  mesh  generation. 
CoBi  can  handle  a  broad  variety  of  grids  including:  simple  Cartesian  grids,  curvilinear 
hexahedra,  Octree  grids,  unstructured  tetrahedral  grids,  and  general  polyhedral  grids.  Unlike  in 
conventional  engineering  where  the  geometry  is  “man  made”  and  defined  by  mathematical 
formulas  (lines,  surfaces,  solids)  cell  biology  is  not  well  defined.  It  is  very  challenging  to 
describe  cell  geometry  in  a  mathematical  form.  The  best  way  is  probably  to  use  a  microscopy 
image  of  a  cell  and  image  processing  software  tools  to  generate  the  geometry  and  to  generate 
mesh  for  that  cell.  Figure  3.1  presents  an  example  approach  developed  at  CFDRC  to  generate  a 
Quadtree  (in  2D)  and  Octree  (in  3D)  mesh  on  a  biological  cells. 


Figure  3.1:  Image  based  generated  cell  geometry  and  a  quadtree  mesh  for  a  bacterial  cell. 

3.3  Membrane  Model 

In  cell  biology,  membranes  play  an  extraordinary  role  protecting  cells  or  organelles  from  the 
outside  environment,  controlling  cell  shape,  enabling  in-transport  of  nutrients  and  excretion  of 
waste  products.  Membrane  bound  receptors  (proteins  and  peptides)  play  the  sensing, 
communications,  and  gate  keeping  roles.  Despite  the  differences  in  membrane  composition  for 
different  cell  types  we  will  treat  them  as  very  thin  shells  with  embedded  receptors,  modeled  as 
either  bound  or  diffusible  “species”.  The  receptors  can  undergo  chemical  reactions  with 
external,  internal,  and  other  membrane  bound  species.  3D  transient  (4D)  surface  diffusion- 
reaction  equations  on  the  membrane  surface  will  be  solved. 

The  membrane  geometry  can  be  time  dependent,  can  move,  and  can  undergo  deformations.  In 
the  present  version,  local  velocities  (deformations)  of  the  membranes  will  be  prescribed  or 
computed  from  other  physiological  models.  In  the  computational  model,  it  is  assumed  that  the 
membrane  surface  can  be  defined  as  a  set  of  interconnected  triangles  forming  an  “unstructured” 
3D  mesh. 

3.3.1  Basic  Equations  and  Assumptions 

Based  on  experience,  the  governing  equations  for  membrane  species  are  described  in  the  integral 
form.  It  ensures  strong  conservation,  easy  formulation  of  moving/deforming  geometry,  and  it 
provides  intuitive  framework  for  numerical  approximation  using  Finite  Volume  Method  (FVM) 
on  arbitrary  control  volume  shapes.  The  time-dependent  scalar  conservation  law  will  be  used  in 
the  following  form: 
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(3.8) 


dt 


oD^nV(p.da  =  |  S^dV 


a 


V 


where  (pi  is  the  i-specie  concentration,  V  is  the  control  volume  (function  of  time),  o  is  the  surface 
of  the  control  volume,  v,  and  Vg  are  the  advection  velocity  and  grid  velocity  (e.g.  due  to  the  mesh 
deformation)  vectors,  n  is  the  surface  normal  vector  for  control  volume  faces,  Di  is  the  diffusion 
coefficient  for  specie  i,  and  Si  is  the  source  term  e.g.  due  to  reaction  rate.  The  equation  is  solved 
in  the  inertial  frame,  so  if  the  entire  body  of  the  membrane  is  moving  v=0. 

For  moving/deforming  grids,  the  volume  conservation  principle  has  to  be  ensured: 

f  J  dV  =  ^v^-Hda  (3.9) 

V  G 

This  formulation  is  general  enough  to  enable  future  extensions  of  the  membrane  model.  Note 
that  to  represent  surface  bound  (immobile  or  non-diffusive)  receptors  it  is  enough  to  set  D=0. 

3.3.2  Numerical  Intesration  Method 

To  numerically  solve  the  governing  equations  FVM  is  used  and  integrate  over  an  arbitrary 
control  volume.  Note  that  the  surface  integrals  are  replaced  by  flux  terms  across  control  volume 
faces.  All  dependent  variables  are  stored  at  the  cell  centers,  mesh  coordinates  at  the  mesh 
vertices,  and  all  fluxes  are  calculated  at  the  cell  faces. 

The  transient  term  is  integrated  with  a  backward  Euler  method: 


d 

dt 


v<pi-v>: 

At 


(3.10) 


The  source  term  is  a  simple  product  of  reaction  rate  times  the  cell  volume.  Note  that  the  reaction 
rate  is  a  nonlinear  function  of  concentrations  (pi  and  which  needs  to  be  linearized.  In  the  current 
CoBi  version,  simple  bimolecular  as  wells  as  Arrhenius,  Michelis-Menten,  and  other  more 
general  expressions  are  used. 

The  diffusive  fluxes  are  being  calculated  as  surface  integrals.  Figure  3.2  illustrates  the  control 
volume  notation  used  for  flux  evaluation. 
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C2 


Figure  3.2:  Nomenclature  used  for  Calculation  of  Cell  Face  Fluxes,  and,  Unstructured  control 

volume  notation. 


For  convenience,  the  details  of  flux  evaluation  of  a  seleeted  faee  denoted  as  e  are  presented.  The 
same  rules  apply  to  all  faees  of  the  eontrol  volume.  The  diffusive  flux  of  (p  is  eomprised  of  two 
terms,  main  and  cross  terms  as  follows: 

De  =  I D^  f2(pn,dF  «  D:  +  (3.11) 


where  Ae  is  the  area  of  face  e  of  the  control  volume,  and  D“  and  represent  main  and  cross 
diffusion  fluxes  through  face  e.  By  discretizing  the  gradient  of  (p  at  the  cell  face  eenter  e  with  the 
central  differenee  scheme,  the  following  expressions  are  obtained: 


07=0 


(p,e 


-{(pE-(Pp) 


(3.12) 


Afn^-iif 


(p,e 


f2-U\ 


^i(Pc2-(Pcx) 


(3.13) 


where  (pci  and  (pc2  are  the  eoncentration  values  at  the  corners.  The  values  at  the  comers  are 
evaluated  from  the  cell  center  values  via  interpolation  (e.g.  1/r  averaged). 

The  adveetive  flux  is  evaluated  as  follows: 


Fe=\(p7-n,dr«{v,^A)^-(p^ 


(3.14) 


where  Vn  is  the  faee  normal  relative  (transport-grid)  veloeity  and  (pe  is  the  eoncentration  vale  at 
the  center  of  the  e  faee.  For  praetical  reasons  we  ealculate  cell  faee  value  based  un  the  “upwind” 
value  from  the  neighbor  cell  center  i.e.: 
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(3.15) 


(Pe 


(Pe---v„<0 


All  terms  are  eomputed  fully  implicitly.  To  minimize  the  computer  storage  requirements  the 
cross  diffusive  fluxes  are  calculated  semi-implicitly  (previous  iteration  value). 

The  grid  vertex  velocity  is  specified  from  another  (LBL  main  solver)  software  module  but  the 
cell  face  velocities  are  computed  within  the  membrane  module  in  CoBi  based  on  the  volumetric 
conservation  law  (shown  above)  as  from  either  Backward  Euler  or  Crank  Nicholson 
approximation  as  follows: 


V 


g 


At 


2At  ^ 


(3.16) 


where  Vg  is  the  grid  velocity  normal  to  the  face,  Vg°  is  the  old  time  step  value,  AVf  is  the  volume 
swept  by  the  moving  face  of  the  deforming  control  volume,  and  At  is  the  time  step. 

After  summing  up  all  cell  face  fluxes  and  volume  integral  terms,  the  discretized  equation  can  be 
expressed  in  the  following  algebraic  form: 

ttpCPp  =  (3-17) 

nb 


Where  a?  is  the  main  diagonal  coefficient  and  anb  are  the  neighbor  link  coefficients,  which 
involve  contributions  of  advective  and  diffusive  fluxes.  This  system  of  algebraic  equations 
forms  a  matrix  with  a  constant  bandwidth  of  only  four  nonzero  entries  for  the  triangle  cells.  One 
of  available  sparse  matrix  solvers  (e.g.  Conjugate  Gradient  Squared,  CGS,  or  GMRES)  is  used  to 
solve  the  equation  set.  It  is  solved  sequentially  for  each  of  the  concentrations. 

3.3.3  Model  of  Membrane  Geometry 

The  membrane  surface  is  represented  in  the  form  of  a  graph  of  triangles,  with  coordinates 
defined  at  the  vertices.  Figure  3.3  illustrates  example  cell  membrane  triangulations.  Each  vertex 
also  has  the  value  of  the  cell  membrane  thickness.  The  flexibility  to  name  surface  patches  is  also 
provided  to  enable  different  volume  and  boundary  condition  specifications,  and  data  post 
processing. 


Figure  3.3:  Cell  membrane  triangulations. 
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In  the  first  CoBi  release,  uniform  thiekness  is  assumed.  The  membrane  module  automatieally 
caleulates  all  geometrieal  parameters  of  for  each  volume  including,  areas,  volume,  surface 
normals,  distances  etc. 

3.4  Biochemistry  Model 

The  source  term  S,  (also  marked  as  R,)  in  the  species  transport  equations  (eq.  3.7)  is  computed 
based  on  the  mechanism  of  reaction  kinetics.  The  net  source  of  chemical  species  i  due  to 
reaction  R  ,  is  computed  as  the  sum  of  the  reaction  sources  over  the  Nr  reactions  that  the  species 
participate  in: 


r=\ 


(3.18) 


where  Mv,  ,  is  the  molecular  weight  of  species  i  and  Ri,r  is  the  rate  of  reaction  (creation  or 
consumption)  of  species  i  in  reaction  r.  The  reaction  mechanism  can  be  written  in  the  form  of 
reversible  reaction  mass  action  law: 

'Zy'u[C,]  |;V,,,[C,]  r=l,NR  (3.19) 

1=1  k'’r  1  =  1 

where  N  is  the  number  of  species  in  the  complete  reaction  mechanism  and  summation  is  over  all 
species  appearing  in  that  reaction,  v^  and  v*^  are  the  stoichiometric  coefficients  for  reactants  and 
products  respectively  and  square  brackets  denote  symbols  of  the  species,  k^  and  k'’  denote 
forward  and  backward  (reverse)  reaction  rates  constants  for  a  giver  reaction  r.  This  formulation 
is  valid  for  reversible  or  irreversible  (k'’=0)  reactions. 

In  CoBi,  reaction  may  occur  in  the  volume  or  on  the  surface  e.g.  cell  membrane. 

The  individual  reaction  step  rates  for  j-reaction  step  coj  is 

N  ■  N  ■ 

’  species  r  species 

n  c;'’  n  C  (3.20) 

i=\  i=\ 

where  k^j  is  the  forward  rate  constant  and  k'(j  is  the  backward  (reverse)  rate  constant  and  Ci  are 
the  species  concentrations  (molar). 

The  molar  rate  of  production  (R>0  or  destruction  R<0)  of  species  i  in  equation  3.18  can  now  be 
written  as: 


Nr  ,  . 

j^producion  ^  £(V,;  -  V0-(Oj 


J  =  i 


(3.21) 
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where  Nr  is  the  number  of  reactions  species  in  reaction  r,  Cj,r  are  the  molar  concentrations  of 
species  (reactants  and  products  respectively),  and  the  symbol  If  denotes  multiplication.  The 
forward  and  reverse  reaction  kinetics  rates  are  supplied  by  the  user.  They  can  be  functions  of 
different  varaibles  (e.g.  temperature).  In  some  cases  the  forward  rate  and  equilibrium  constant 
Kr  are  known  allowing  the  reverse  rate  to  be  simply  computed  as: 


k\  = 


(3.22) 


Another  type  of  biochemical  reactions  in  CoBi  are  the  enzyme-catalyzed  reactions  in  which  the 
enzyme  E,  a  catalyst,  and  substrate  S  participate  in  a  reversible  reaction  forming  enzyme- 
substrate  complex  ES  followed  by  an  irreversible  reaction  of  forming  the  product  P  and  recycling 
the  enzyme.  The  reaction  mechanism  can  be  expressed  as: 

k^ 

S  +  E<^^ES^E  +  P  (3.23) 

k,i 


The  above  enzyme  catalyzed  kinetics  reaction  system  can  be  expressed  as  either  multi  step  mass 
action  kinetics  or  Michaelis-Menten  kinetics.  For  enzymatic  reactions  in  closed  systems,  kinetic 
scheme  (3.23)  implies  the  following  conservation  relations: 


Et=E  +  ES,  (3.24) 

and 

St=  S  +  ES  (3.25) 

where  Ej-  and  Sj-  are,  respectively,  the  total  concentrations  of  the  enzyme  and  the  substrate  in  the 
system.  Applying  the  law  of  mass  action  to  kinetic  scheme  (3.23),  and  using  equation  (3.24- 
3.25)  to  eliminate,  yields  the  following  set  of  differential  equations: 

Try 

—  =  -k^-(E^-ES)-S  +  k^-ES  (3.26) 

dt 


dES 

dt 


=  k,-{{E^-ES)-S-K^-ES 


(3.27) 


—  =  k^-ES 
dt  " 


where  the  Michaelis  constant  Km  is  defined  as 


Km  = 


k^^  -h  k2 

K 


(3.28) 


(3.29) 
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This  reaction  scheme  (3.23)  is  mathematically  described  by  a  set  of  coupled  non-linear  second- 
order  algebraic  (3.24-3.25)  and  differential  equations  (3.26-3.29)  that  could  be  numerically  stiff. 
Stiff  ODE  solvers  are  then  needed.  When  a  quasi-steady-state  assumption  is  made  the  initial  rate 
of  product  formation  can  be  expressed  as  a  function  of  initial  substrate  concentration,  and  then 
obtaining  kinetic  parameters  by  solving  the  Michaelis-Menten  (MM)  equation: 


dt  " 


(3.30) 


where  v^ax  is  the  maximum  velocity  (product  formation  rate)  and  So  initial  substrate 
concentration. 


Michaelis  Menten  law  is  also  known  as  saturable  or  capacity-limited  reactions  and  can  be 
simplified  for  limiting  conditions  for  the  substrate  concentration  S.  If  S»Km  then  it  reduces  to 
Vmax(  i.e.  zero  order  reaction).  The  enzyme  catalyzed  kinetics  reaction  system,  Eq.  3.23,  has 
been  extended  to  more  complex  biochemical  reaction  systems  involving  multiple  substrates, 
multiple  products,  enzyme  inhibitors  and  activators.  Table  3.1  presents  examples  of  such 
reactions. 


Table  3.1:  Examples  of  enzyme  catalyzed  kinetics  reactions  and  reaction  rate  expressions,  Cs 
and  C]  are  the  substrate  and  enzyme  concentrations  while  co  is  the  reaction  rate. 


Type 

Simplified  Kinetics 

Stepwise  Rate  Expression 

Michaelis-Menten  - 
MM 

Km 

S^P 

V 

^  max 

V  c 

00=  ^ 

Cs+K 

MM  plus  inhibition 
and  activation 

K 

S  ^  P 

ke^ 

^  .  Kx  Cs 

1  +  —  -1-  — 

Q  K, 

Competitive 

Inhibition 

keo 

S  +  I  ^  P  +  I 

Ks,K, 

ke^Co 

ao  = 

I  Cl 

Cs+Ks  1  +  ^ 

V  J 

Non-Competitive 

Inhibition 

kCo 

S  +  I  ^  P  +  I 

Ks,K, 

ke  Cs 

0)  =  f  \ 

(Q  +  1  + 

V  ) 

Multiple 

Substrate/One 

Enzyme 

ki,k2,eo 

S,  +  Sy  ^2 

^1,^2 

1  Ke^Cs^ 

CO  =  - 7 - i - : 

r  Cs  Cs 
K,  \+  ^‘  + 

1  K, 

\ 

/ 
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In  several  cases  one  needs  to  model  ligand  L  and  receptor  R  binding  e.g.  on  cell  membranes.  A 
simple  L-R  binding  can  be  described  as  a  reversible  reaction: 


R  +  L< 


k 


r 


(3.31) 


Where,  R  is  the  membrane  receptor  available  to  bind  with  free  ligand  (e.g.  drug)  L  and  C  is  the 
bound  ligand-receptor  complex  drug,  and  kf,  and  kr  are  forward  and  backward  rate  constants. 
The  rate  of  change  in  concentration  of  bound  drug  and  fee  drug  is  defined  as: 


—  =  kf-R-L-kC 
dt  ^ 

(3.32a) 

—  =  -kf-R-L  +  kC 
dt  ^ 

(3.32b) 

—  =  -kf-R-L  +  kC 
dt  ^ 

(3.32c) 

Assuming  kf  and  kr  are  known,  the  three  equations  for  R,  L,  and  C  can  be  solved.  Assuming  that 
there  is  a  finite  number  of  receptor  binding  sites  and  introducing  the  maximum  cell  binding 
receptors,  Bmax,c  so  that 


5.ax=^  +  C  (3.33) 

Substituting  equation  (3.33)  into  (3.32a)  to  remove  the  receptor  R,  there  are  only  two  equations 
for  free  and  bound  drug.  For  example,  the  equation  for  L  is: 

^  =  (3.34) 

The  biochemistry  module  in  CoBi  solves  systems  of  Ordinary  Differential  Equations  (ODEs) 
individually  for  each  control  volume  in  the  3D  and  membrane  surface  computational  domain. 
The  solution  of  ODE  equations  provides  species  concentrations  for  the  PDE  solver  module. 
They  are  used  to  compute  reaction  kinetics  source  terms  in  the  PDE  equations. 

3.5  Boundary,  Initial,  and  Volume  Conditions 

The  mesh  imported  to  CoBi  has  to  have  following  data: 
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•  Number  of  nodes,  faces,  and  volumes, 

•  X,Y,Z  coordinates  for  all  nodes, 

•  Node  link  lists  to  cells, 

•  Boundary  condition  path  names  and  face  link  lists, 

•  Volume  condition  zone  names  and  cell  lists. 

The  names  (tags)  of  boundary  and  volume  conditions  are  used  to  impose  Boundary  Conditions 
(BCs)  and  Volume  Conditions  (VCs)  on  these  surface/volume  meshes.  BCs  specify  the 
dependent  variables  at  the  domain  boundaries  whereas  VCs  specify  material  types  and  models  to 
be  used  in  different  parts  of  the  computational  domain.  At  present  CoBi  users  can  specify  them 
in  the  *.SIM  fde  (a  text  file  used  as  input  to  run  CoBi)  or  interactively  via  JCoBi  GUI  (as 
described  in  JCoBi  User’s  Manual).  Note  that  JCoBi  creates  the  *.SIM  file  so  the  user  can 
verify  or  modify  the  BC  and  VC  types  and  property  values  using  any  text  editor.  There  may  be 
different  VCs  for  cell  cytoplasm,  for  extracellular  space  and  for  the  membrane. 

The  boundary  types  available  in  CoBi  are  as  follows: 

•  Fluid  flow  inlet  and  exit:  with  specified  pressure,  velocity,  or  mass  flow 
.  Wall 

•  Symmetry 

•  Periodic  boundary, 

•  Internal  faces  e.g.  membranes 

3.6  Material  Properties 

CoBi  simulation  requires  specification  of  models  (e.g.  flow,  mixing,  reaction,  etc.)  and  the  setup 
of  the  physical  properties  of  the  material.  The  user  needs  to  input  all  material  properties. 
Material  properties  are  associated  with  volume  conditions  (VCs)  which  are  “tagged”  by  names  in 
the  grid  file.  CoBi  user  inserts  the  values  of  specific  properties  in  the  *.SIM  file  (a  text  file  used 
as  input  to  run  CoBi)  or  interactively  via  JCoBi  GUI  (as  described  in  JCoBi  User’s  Manual). 
Note  that  JCoBi  creates  the  *.SIM  file  so  the  user  can  verify  or  modify  the  BC  and  VC  types  and 
property  values  using  any  text  editor 

Typical  properties  needed  for  CoBi  simulations  are: 

•  Density  and/or  molecular  weights 

•  Viscosity 

•  Diffusion  coefficients 

Properties  may  be  temperature  and/or  composition-dependent.  The  material  properties  are  setup 
based  on  volume  condition  (VC)  name  tags  associated  with  the  grid. 
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4.0  DESCRIPTION  OF  GOBI  CODE  FRAMEWORK 


4.1  Overall  Software  Framework 

CFDRC  has  developed  core  solver  modules  and  the  integrated  software  framework  for 
computational  biology.  It  has  evolved  from  an  initial  membrane  only  solver  to  the  coupled 
volume-membrane  convection-reaction-diffusion  equations  coupled  to  the  ODE  biochemistry 
module.  The  integrated  solver  module,  CoBi,  has  been  tested  on  several  problems  related  to  cell 
biology  and  was  delivered  to  DARPA  (SRI)  in  February  2005.  This  section  briefly  presents  the 
technical  capabilities  and  software  structure  of  CoBi  tools. 

The  CoBi  solver  can  solve  flow  and  diffusion-reaction  problems  in  any  geometry  ranging  from 
single  OD  compartment,  to  ID  plug  reactor,  to  network  of  compartments  connected  by  channels, 
to  ID,  2D,  or  3D  geometries  of  a  single  cell  or  multi-cellular  tissues.  The  unique  capability  of 
CoBi  is  that  in  addition  to  volumetric  equations,  it  can  simultaneously  solve  PDE  reaction- 
diffusion  equations  on  thin  membranes  (surfaces)  surrounding  biological  cells  or  subcellular 
organelles,  endothelial  or  epithelial  walls  of  vessels.  The  geometry  may  change  in  time,  may 
deform,  as  part  of  the  solution  process.  For  example,  the  cell  can  grow  during  the  cell  division 
process.  The  code  has  the  following  modeling  capabilities: 

•  Steady-state  or  transient  analysis 

•  Incompressible  Newtonian  flows 

•  Heat  transfer, 

•  Chemical  species  mixing  and  reaction  in  volume  and  on  membranes 

•  Arbitrary  volumetric  sources  of  mass,  momentum,  and  chemical  species 

•  Flow  and  species  diffusion-reaction  in  porous  media 

CoBi  solver  internally  uses  unstructured  meshes  but  can  impart  structured,  multilink  structured, 
curvilinear  body  fitted,  tetrahedral,  octree,  and  any  finite  element  mesh. 

4.1.1  Preprocessins,  Geometry,  Meshins,  and  Model  Setup 
At  present  CoBi  reads  two  main  files: 

a)  Mesh  data  with  named  volumetric  regions  that  are  used  to  set  boundary  and  initial 
conditions  as  well  as  volumetric  properties  (diffusivity,  density,  porosity,  etc.),  and 

b)  Text  file,  *.SIM,  providing  problem  definition  specification  including  values  for  boundary 
conditions,  material  properties,  numerical  control  parameters,  etc. 

A  prototype  simple  geometric  modeler  and  mesh  generator  have  also  been  developed,  and 
implemented  it  in  the  CoBi  GUI  module,  JCoBi,  for  problem  setup.  To  enable  an  efficient  setup 
of  biochemistry  equations,  direct  interfaces  to  the  open  source  cell  biochemistry  modeling 
languages.  System  Biology  Modeling  Language  (SBML),  was  developed.  It  allows  the  user  to 
import  cell  biology  models  setup  with  other  Bio-SPICE  tools  such  as  IDesigner,  Jamac,  and 
Jigcell. 
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4.1.2  CoBi  Software  Pro2rammin2 


A  general  and  modular  C++  software  framework,  CoBi,  for  solving  advection-diffusion-reaction 
partial  differential  equations  (PDEs)  on  unstructured  grids  was  developed.  The  CoBi  software 
has  been  written  in  C++  using  object-oriented  programming  with  very  powerful  database  for 
handling  geometry,  mesh,  solution  fields,  properties,  and  boundary/volume  conditions.  It 
combines  the  full  Navier  Stokes  equation  flow  solver  for  incompressible  flows,  and  a 
comprehensive  model  of  volume  and  membrane  surface  constrained  convection-diffusion- 
reaction  model.  It  has  database  interfaces  to  geometry  meshing  tools  and  a  direct  interface  to  3D 
visualization/animation  tool  Visit  (www.llnl.gov/visit/).  Figure  4.1  presents  the  complete  CoBi 
toolkit  environment  for  modeling  cell  and  tissue  biology. 


CoBi  Tools  Simulation  Environment 


JigCell 

Tools 


>ceilML 

Systems  Biology 
Markup  Language 


Figure  4.1:  Cell  biology  simulation  environment  and  the  CoBi  modeling  toolkit. 

4.1.3  Visualization  with  Visit  Tools 

Visit  is  a  free  (http://www.llnl.gov/visit/about.htmlj  interactive  parallel  visualization  and 
graphical  analysis  tool  for  viewing  scientific  data  on  UNIX  and  PC  platforms.  Visit  was 
developed  by  the  Lawrence  Livermore  National  Laboratory  (LLNL),  Department  of  Energy 
(DOE)  Advanced  Simulation  and  Computing  Initiative  (ASCI)  to  visualize  and  analyze  the 
results  of  terascale  simulations.  Users  can  quickly  generate  visualizations  from  their  data, 
animate  them  through  time,  manipulate  them,  and  save  the  resulting  images  for  presentations.  It 
can  be  used  to  visualize  scalar  and  vector  fields  defined  on  two-  and  three-dimensional  (2D  and 
3D)  structured  and  unstructured  meshes.  Visit  was  designed  to  handle  very  large  data  set  sizes 
in  the  terascale  range  and  yet  can  also  handle  small  data  sets  in  the  kilobyte  range. 
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Visit  development  at  LLNL  began  in  the  summer  of  2000,  and  the  initial  version  of  Visit  was 
released  in  the  fall  of  2002.  Visit  supports  C++,  Python  and  Java  interfaces.  The  C++  and  Java 
interfaces  make  it  possible  to  provide  alternate  user  interfaces  for  Visit  or  allow  existing  C++  or 
Java  applications  to  add  visualization  support.  This  feature  will  allow  CFDRC  to  adapt  VisIT  as 
a  Cell  Biology  Visualizer.  Figure  4.2  presents  an  example  GUI  window  from  the  VisIT 
package. 


Visit  Main  window 


_Vis  window 


Help]  lB  g^UlP'*”  LU  m  yiii  I  III  ^1  Xr°  ► 


Pseudocolor  Plot 
Vo^pro^ljyjc 


Figure  4.2:  GUI  Window  of  the  VisIT  visualization  toolkit  used  for  modeling  cell  biology  in 

CoBi. 

4.2  CoBi  Solver  Structure 

The  CoBi  solver  is  written  in  model  object  oriented  C++  software  compiled  on  Linux  and 
Windows.  It  is  written  in  the  form  to  enable  future  addition  of  new  modules  for  physics  and 
biology  disciplines  such  as  electrokinetics  for  neuronal  cells,  stress/deformation  for  modeling 
cell  biomechanics,  and  other  modules.  The  CoBi  code  consists  of  the  following  modules: 

In-core  database  with  all  variables,  properties,  mesh,  solution 
I/O  routines  to  input  files,  to  mesh,  and  to  visualization  tools 
Main  driver 

Mesh  and  geometry  initialization  module 

Physical  data  initialization 

Mesh  movement  module 

Flux  calculation  module  (diffusive,  advective) 

Flow  module 

Reaction  diffusion  module 
Thermal  module 

Simple  chemical  kinetics  source  terms  evaluation  module  (including  linearization) 

Matrix  assembly 


27 


Linear  equation  solvers 
CoBi  software  framework  eonsists  of 

•  JCoBi  a  preproeessor  graphical  user’s  interface  (GUI)  for  interactive  setup  of  problem 
biophysics,  boundary  and  volume  conditionings  setup,  specification  of  material 
properties  and  solution  control  parameters.  An  interface  to  biochemical  pathway  tools 
such  as  JDesigner  is  enabled  via  the  SBML  interface, 

•  CoBi  solver,  the  main  computational  engine  solving  all  the  equations,  and 

•  Visit  -  a  graphical  post  processing  for  3D  results  visualization  and  animation. 

At  present,  the  geometry  modeling  and  grid  generation  module  has  limited  capabilities.  JCoBi 
can  setup  2D  and  3D  geometries  for  typical  cell  shapes  such  as  sphere,  ellipsoid,  box,  etc. 
However,  the  code  allows  the  user  to  import  an  externally  generated  grid  fde  with  annotated 
boundary  and  volume  conditions.  One  potential  approach  to  generate  2D  cell  biology  meshes  is 
to  use  the  CFDRC  Micromesh  tool.  The  next  section  describes  this  approach. 

JCoBI  is  written  in  Java  while  the  CoBi  solver  is  written  in  the  C++  language.  CoBi  tools  can  be 
used  on  computers  with  either  Windows  or  Linux  operation  systems. 

Figure  4.3  shows  the  organizational  structure  of  CoBi  modeling  framework. 


Figure  4.3:  Schematics  of  CoBi  Program  Structure. 

4.3  Automated  Geometry  and  Meshing  of  Evolving  Biological  Cells 

During  this  project,  CFDRC  has  adapted  the  CFD-Micromesh  software  tool  for  the  application  in 
automated  image  processing,  geometry  construction,  and  mesh  generation  for  2D  and  quasi  3D 
simulation  of  cellular  biophysics  and  biochemistry.  Figure  4.4  presents  example  results  in  terms 
of  original  image,  unstructured  triangular  mesh,  adaptive  binary  tree  mesh,  and  multi  block 
structure  mesh  for  the  Bacillus  Subtilis  cell  during  germination  process.  Note  the  cell  membrane 
capturing  capability  and  formation  of  septum. 
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Figure  4.4\  Bacillus  Subtilis  during  germination  cell  cycle  division.  Image  and  mesh  examples 

obtained  with  CFD-Micromesh. 


4.4  Preprocessing  and  Geometry/Mesh  Generation 


At  present,  the  geometrical  model  preprocessing  is  achieved  by  either  using  JCoBi  for  pre  set 
parametric  geometries  of  selected  cell  shapes,  or  reading  an  external  text  fde  (*.SIM)  for  the 
model  setup  and  a  DTP  fde  for  the  computational  mesh.  In  the  second  case,  the  mesh  has  to  be 
generated  by  other  software  tools.  JCoBi  provides  simple  but  complete  capability  to  generate 
typical  geometries  representing  cell  shapes  such  as  a  sphere,  rod,  ellipsoid,  cube,  and  their 
parametric  variations.  For  block  shapes,  the  user  can  selected  structured  mesh  while  the 
spherical  and  rod  shape  geometers  are  meshed  with  an  unstructured  3D  tetrahedral  mesh.  The 
Figure  4.5  presents  a  Geometry/Mesh  Generation  window  in  JCoBi. 


Figure  4.5:  JCoBi  GUI  Interface  for  generating  3D  cell  biology  geometry. 

JCoBi  provides  graphical  capability  to  set  boundary  and  volume  conditions.  An  imported  mesh 
should  have  name  tags  attached  to  boundary  condition  mesh  patches  and  volume  conditions 
name  tags  to  specific  volumes.  By  “clicking”  a  boundary  condition  name  tag  in  the  GUI 
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window  JCoBI  will  highlight  that  boundary  and  user  can  set  the  boundary  condition  values. 
Similar  capability  is  available  for  setting  volume  condition  values  such  as  material  properties 
and  for  activation  of  physical  models.  Figure  4.6  presents  a  JCoBi  GUI  snapshot  from  a 
window  for  setting  boundary  and  volume  conditions. 


Figure  4.6\  JCoBi  GUI  window  for  setting  boundary  and  volume  conditions. 

To  simulate  a  cell  biology  problem  with  CoBi  the  user  needs  to  setup  a  “simulation  fde” 
containing  information  about  mesh,  boundary  conditions,  material  properties,  numerical  solution 
controls,  etc.  The  model  setup  fde  (fdename.SIM)  may  be  generated  manually  through  a  text 
editor,  or  generated  through  the  newly  developed  JCoBi  GUI  module,  for  a  user  friendly 
graphical  problem  setup.  Also  from  the  GUI,  the  simulation  may  be  run  and  monitored  for 
convergence  as  described  in  the  next  section. 

4.5  Post-processing 

The  computational  results  obtained  with  CoBi  can  be  stored  in  two  formats: 

DTF  developed  by  CFDRC  and  freely  available  form  CFDRC  web  site,  and 
VTK  an  open  source  database  format  used  for  imaging. 

For  advanced  3D  visualization  and  animation  of  simulation  results  we  recommend  using  Visit 
tools  for  the  visualization  of  computational  results.  For  more  details  see  Section  4.1.3. 

For  the  preliminary  analysis  of  simulation  results  the  user  can  use  JCoBi  tools.  JCoBi  provides 
capabilities  to  see  3D  surfaces  with  color  mapped  solution  variables  e.g.  specie  concentration  or 
temperature.  During  the  execution  of  a  CoBi  ’’job”  the  user  can  monitor  the  progress  of  the  job 
using  the  JCoBi  line  plotting  capability.  It  allows  the  user  to  display  (plot)  the  convergence 
characteristics  and  values  of  dependent  variables  at  a  monitoring  point.  Figure  4.7  presents 
example  JCoBi  windows  with  a  3D  surface  mapping  of  a  specie  concentration  of  bacterial  rod 
geometry  and  an  x-y  plot  of  a  time  course  of  a  dependent  variable  (specie  concentration)  at  a 
monitoring  point.  The  detailed  capability  of  the  JCoBi  tool  can  be  found  in  the  CoBi  User’s 
Manual. 
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Figure  4. 7:  JCoBi  GUI  interface  for  3D  graphical  post  processing  of  a  species  concentration  of 
a  bacterial  surface  and  an  x-y  plot  of  a  time  course  of  a  dependent  variable  (specie 

concentration)  at  a  monitoring  point. 

4.6  Example  Problems  Setup  in  JCoBi  and  Execution  in  CoBi 

This  section  presents  a  brief  description  of  the  simulation  process  in  the  CoBi  framework.  An 
example  problem  of  a  three  dimensional  simulation  of  blood  flow  through  a  channel  with  two 
square  cells  located  in  the  flow  field,  forcing  flow  around  the  two  cells  is  demonstrated. 

GEOMETRY :  A  three  dimensional  geometry  and  grid  is  generated  for  the  simulation.  After 
importing  the  mesh  into  JCoBi  the  user  can  setup  the  associated  volume  and  boundary  conditions 
as  shown  in  Figure  4.8. 


Outlet  Boundary 
“Blood  Outlet” 


“Blood” 

Volume 


“Cell2”  Volume 
and  “Blood_Cell2” 
Boundary 


Boundary 
“Blood  Inlet” 


“Celll”  Volume  and 
“BloodCelll”  Boundary 


Top,  Bottom  and 
Side  Boundaries 
“AIR  Blood” 


Figure  4.8:  JCoBi  GUI  interface  for  setting  boundary  and  volume  conditions. 
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PROBLEM  SETUP: 


JCoBi,  the  graphical  user  interface  for  CoBi, 
may  be  used  to  setup  the  simulation.  The 
main  window  upon  opening  the  JCoBi  is 
shown  in  Figure  4.9.  If  a  simulation  file 
(*.SIM)  already  exists,  it  may  be  opened 
through  the  “Menu”  pull-down  menu,  or  a 
new  problem  may  be  setup  and  saved. 


Figure  4.9\  Main /CoRz  window. 


Flow  Setup:  _ 

[f  Modules  j'^lunies  Boundaries~Y  Simulation  | 

E  Flow  Module 
□  Fleat  Module 


Under  the  “Modules”  tab  (Figure  4.10), 

activate  the  flow  module  by  checking  the  box  i - 

next  to  the  “Flow  Module”.  I - . 

Edit  Fteactions 
Add  Species 


Figure  4.10:  Module  selection. 


Volume  Conditions: 

The  volumes  presented  in  the  simulation  are 
listed  under  the  “Volumes”  tab  (Figure  4.1 1). 
By  selecting  each  volume  in  the  list,  the 
necessary  settings  may  be  entered  for  each 
volume.  After  completing  the  settings  for 
each  volume,  click  the  “Apply”  button  to 
accept  the  settings.  The  following  table 
summarizes  the  required  settings  for  this 
simulation,  and  the  image  illustrates  the 
settings  for  one  of  the  volumes.  Note  that 
since  there  is  to  be  no  flow  in  the  two  cells,  a 
viscosity  value  of  zero  is  applied  to  these 
volumes. 


f  Modules  \  Volumes  \  Boundaries  Simulation  | 

Volumes  List: 

Blood 

Celll 

Cell2 


Figure  4.11:  Volume  selection. 
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Table  4. 1 :  Volume  conditions. 


Volume 

Density 

Reaction 

Viscosity 

Blood 

1025.0 

None 

0.00115 

Celll 

1125.0 

None 

0 

Cell2 

1125.0 

None 

0 

Boundary  Conditions: 


The  boundary  conditions  presented  in  the 
simulation  are  listed  under  the  “Boundaries”  tab 
(Figure  4.12).  By  selecting  each  boundary  in  the 
list,  the  necessary  settings  may  be  entered  for  each 
boundary.  After  completing  the  settings  for  each 
boundary,  click  the  “Apply”  button  to  accept  the 
settings.  The  following  table  summarizes  the 
required  settings  for  this  simulation,  and  the  image 
illustrates  the  settings  for  one  of  the  boundaries. 
Note  that  for  boundaries  where  there  is  to  be  no 
flow,  velocity  values  of  zero  are  used  (u,  v,  w 
velocity  components  are  required). 


( Modules  Volumes  Boundaries  Simulation  | 

Boundaries  List: _ 

AIR_Blood _ 

^loodjnlet 

BlDDd_Outlet 

Blood_Cell1 

Blood_Cell2 


Figure  4.12:  Boundary  selection. 


Table  4.2:  Boundary  conditions. 


Boundary 

BC  Type  for  Flow 

Value 

AIR  Blood 

Velocity 

0  0  0 

Blood  Inlet 

Velocity 

0  l.OE-4  0 

BloodOutlet 

Pressure 

0 

Blood  Celll 

Velocity 

0  0  0 

Blood  Cell2 

Velocity 

0  0  0 

Simulation  Settings: 


The  simulation  control  settings  can  be  set  under  the 
“Simulation”  tab  (Figure  4.13).  The  simulation 
control  settings  for  this  simulation  are  illustrated 
below.  By  default,  DTF  output  files  will  be 
generated  unless  VTK  output  files  are  selected. 


f  Moilules~~(  Volumes  \  Boundaries  (  Simulatinn 

Iteration: 

Number  of  steps: 

Time  step  [s]: 

Dump  interval: 

Dump  to  VTK  files: 

Moving  Grid: 


lio 


25 

0.05 

1 _ 

□ 

□ 


Run  Simulation 


Show  Plots 


Figure  4.13:  Simulation  control  settings. 


33 


Run  Simulation: 


Before  running  the  simulation,  save  the  simulation  under  the  pull-down  “Menu”  button  in  JCoBi. 
Once  the  file  is  saved,  click  the  “Run  Simulation”  button  under  the  “Simulation”  tab  in  JCoBi. 
Convergence  of  the  simulation  may  be  monitored  by  clicking  the  “Show  Plots”  button. 

RESULTS: 

The  simulation  results  generated  by  CoBi  are  stored  in  either  DTE  or  VTK  format.  The  user  can 
import  them  to  any  graphical  post  processing  tools  such  as  Visit.  Figure  4.14  presents  example 
flow  distribution  from  this  demonstration  simulation. 


Figure  4.14:  CoBi  Simulation  results  displayed  in  an  independent  graphics  visualization  tool 

such  as  VisIT. 
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5.0  RESULTS  OF  CELL/TISSUE  BIOLOGY  SIMULATIONS  WITH  CoBi 


In  CFDRC's  previous  BioCOMP  Program,  a  subcontract  project  with  LBL,  CFDRC  has 
developed  stand  alone  software  for  spatio-temporal  modeling  of  biological  membranes.  A  3D 
cell  biology  Membrane  Model  has  been  developed  and  delivered  in  source  form  to  LBL  and  SRI. 
Summary  of  that  project  and  example  simulation  results  are  presented  in  the  Appendix  A  to  this 
report.  In  this  project  the  Membrane  Model  has  been  used  as  the  starting  point  for  the 
development  of  the  CoBi  code.  This  section  presents  results  obtained  with  the  CoBi  code  for 
biological  cell  problems  ranging  from  2D  and  3D  membrane  transport,  to  compartmental  and 
spatiotemporal  model  of  cellular  signaling,  to  3D  tissue  biology  including  perfusion  and 
metabolism. 

All  the  test  cases  presented  in  this  section  have  been  obtained  with  the  CoBi  code.  The  input 
data,  in  the  form  of  *.SIM  text  files,  used  to  run  these  cases  along  with  technical  comments  and 
explanations  are  provided  in  Appendix  B  to  this  report.  For  an  illustrative  purpose,  the  input 
data  for  the  first  test  case,  membrane  diffusion  problem,  is  provided  chapter  5.1. 

After  completed  simulations,  CoBi  creates  an  output  file  for  visualization  of  the  solution  results. 
The  output  file  is  created  in  the  vtk  format  (www.vtk.org/pdf/file-formats.pdf ).  The  vtk  file 
can  be  viewed  using  Visit  visualization  tools  developed  by  the  Lawrence  Livermore  National 
Laboratory  (www.llnl.gov/visit/),  or  other  visualization  software. 

5.1  Diffusion  on  a  Planar  Membrane 

Diffusion  Problem  on  Planar  Square  Thin  Film  - 

Several  ID,  2D,  and  3D  species  diffusion  problems  have  been  tested  in  CoBi  to  test  code 
accuracy,  robustness,  and  consistency.  For  example  ID  problems  have  been  setup  as  2D  or  3D 
to  test  solution  uniformity,  or  comparison  between  2D/3D  structured  grid  and  unstructured  grid 
results  have  been  made  to  check  code  mesh  independence.  A  simple  one  dimensional  specie 
diffusion  problem  was  chosen  to  validate  the  accuracy  of  CoBi  membrane  model.  The  geometry, 
grid,  and  problem  definition  are  shown  in  the  Figure  5.1.  It  is  a  square  box  with  unit  length  and 
its  thickness  is  set  to  .01  and  the  mesh  is  an  unstructured  triangular  grid.  The  exact  solution  to 
this  problem  is  a  linear  distribution  of  specie  concentration  from  zero  to  one  in  the  horizontal 
direction.  As  expected  and  shown  in  Figure  5.1  CoBi  membrane  model  gives  an  exact  solution 
in  this  case.  Exactly  the  same  results  were  obtained  with  the  structured  mesh. 
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Geometry/Mesh 


Species  Distribution 


Species  as  spatial  function 


Figure  5.1:  Computational  results  for  a  specie  diffusion  problem  on  planar  square  membrane. 
Input  File 

Two  text  based  input  files  are  needed  to  run  the  model: 

Membrane2D.MFG  and, 

Membrane2D .  SIM . 


The  MFG  is  a  mixed  element  format  used  by  CFDRC.  It  contains  grid,  volume  condition,  and 
boundary  condition  information.  The  detailed  description  of  the  format  is  available  from 
CFDRC.  The  *.SIM  file  is  an  input  text  file  created  by  the  user.  It  can  be  created  by  any  test 
editor  (including  Microsoft  Word)  and  contains  information  about  control  parameters,  values  of 
volume  conditions  and  boundary  conditions.  The  line  numbers,  01,  02,  ...  (shown  below)  not 
needed  in  the  *.SIM  file  read  by  CoBi.  Here  they  are  used  to  explain  individual  commands  of 
the  *.SIM  file.  The  *.SIM  file  for  this  problem  is  as  follows: 

01  thickness  0.01 

02  iteration  10 

03  module  share 

04  VC  Membrane  const_dens  1 

05  module  species  fl 

06  sweeps  500 

07  relaxation  0. 1 

08  VC  Membrane  const_diff  0. 1 

09  be  Left  fix_value  0.0 

10  be  Right  fix  value  1.0 

11  be  Top  fix  flux  0.0 

12  be  Bottom  fix  flux  0.0 

13  be  Membrane_outside  fix_flux  0 

14  be  Membrane  inside  fix  flux  0 
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In  the  input  file  keyword  thickness  is  the  membrane  thickness.  Keyword  iteration  is  the  number 
of  iterations  for  the  steady  state  case  or  each  time  step.  Keyword  module  is  used  to  specify  input 
data  for  the  module.  Keyword  share  means  that  input  for  shared  data,  such  as  density. 

Keyword  vc  is  used  to  specify  a  volume  condition  followed  by  a  volume  name,  constant  density, 
and  the  value.  Keyword  species  is  used  for  input  of  species  module.  Keyword  sweeps  is  used  to 
specify  the  maximum  number  of  sweeps  for  the  species.  Keyword  relaxation  is  used  to  input 
linear  relaxation  between  a  previous  value  and  a  new  value.  The  next  line  vc  sets  the  constant 
diffusivity  for  volume  Membrane.  Keyword  be  means  the  input  for  boundary  condition  followed 
by  boundary  name,  boundary  condition  type  (fixed  value/fixed  flux). 

Line  01  sets  membrane  thickness  to  0.01; 

Line  02  sets  number  of  sub-iteration  to  10; 

Line  03  declares  the  start  of  share  module; 

Line  04  sets  density  to  1  for  membrane; 

Line  05  declares  the  start  of  species  fl  equation; 

Line  06  sets  number  of  iterations  for  linear  equation  solver; 

Line  07  sets  relaxation  factor  for  fl  to  0.1; 

Line  08  sets  the  volume  condition  (vc)  of  species  diffusivity  to  0.1  for  membrane; 

Line  09-10  set  fix  concentration  0  on  the  Left  and  Ion  the  Right  Boundaries; 

Line  1 1-12  set  zero  gradient  concentration  on  bottom/bottom  surfaces; 

Line  13-14  set  no  flux  boundary  condition  to  membrane  inside/outside  surfaces; 

Output  File 

For  visualization,  the  solution  is  output  to  files  in  the  vtk  format. 

5.2  Diffusion  on  a  Spherical  Cell  Membrane 

This  test  problem  is  designed  to  validate  and  demonstrate  the  capability  of  the  CoBi  membrane 
model  to  simulate  the  diffusion  of  a  species  constrained  on  a  curved  surface.  It  is  a  challenging 
problem  of  solving  3D  convection-diffusion-reaction  transport  equation  on  a  general  surface  3D 
mesh.  Figure  5.2-a  shows  the  spherical  surface  unstructured  triangulated  grid  and  three- 
dimensional  species  concentration  distribution.  At  the  boundary  conditions  we  set  the 
concentration  of  the  species  equal  to  1  on  a  circular  patch  on  the  top  of  the  sphere,  while  the 
concentration  at  the  bottom  circular  patch  is  set  to  0.  The  sphere  has  a  radius  of  unity  and  the 
membrane  thickness  is  0.01. 

The  front  view  of  the  sphere  with  species  concentration  contour  map  is  shown  Figure  5.2-b.  One 
can  see  that  the  species  has  high  concentration  gradient  near  the  top  and  bottom  regions  due  to 
smaller  diffusion  areas.  The  contour  lines  are  straight,  which  means  the  numerical  method 
preserves  top-bottom  symmetry  and  horizontal  uniformity  while  the  mesh  is  unsymmetrical. 

The  annotated  CoBi  input  file  for  this  test  case  is  provided  in  Appendix  B. 
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5.3 


a.  Geometry  and  Grid  b.  Species  Distribution 

Figure  5.2:  CoBi  membrane  model  for  specie  diffusion  on  a  thin  spherical  surface. 

Diffusion-Reaction  of  a  Flat  Membrane 


This  test  case  is  designed  to  test  the  transient  interaction  between  diffusion  and  chemical  reaction 
between  two  species  on  a  2D  surface.  The  geometry  of  this  case  is  the  same  as  test  case  1 
presented  in  section  5.1  above.  An  unstructured  mesh  has  been  used.  Two  species,  A  and  B,  are 
solved  in  transient  mode  with  a  simple  one  step  chemical  reaction  as  follows: 

d  -h  B  2B  "1“  B 


Both  A  and  B  species  concentrations  are  initialized  to  zero  on  the  surface  (in  this  case  volume 
condition).  Specie  A  =1  is  set  on  the  right  side  boundary  and  B  =1  is  set  on  left  side  boundary 
and  zero  on  the  opposite  sides  respectively.  Figure  5.3a  shows  time  color  contour  map  of  specie 
A  at  a  t=0.6s.  Figure  5.3b  shows  a  sequence  of  profdes  of  species  concentration  distributions 
along  the  x  coordinate  at  t=  0.5s  to  t=ls.  Note  that  until  t=0.7s  both  species  interdiffuse  and 
there  is  little  reaction  (production  of  specie  B).  At  t=ls  specie  B,  being  produced  from  A, 
dominates.  The  results  match  our  expected  trends.  The  annotated  CoBi  input  fde  for  this  test 
case  is  provided  in  Appendix  B. 


Figure  5.3A:  Color  contour  map  of  specie  A. 
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Figure  5.  SB:  Species  A  and  B  concentration  distributions  as  a  function  of  time. 

5.4  Diffusion-Reaction  on  a  Bacterial  Membrane 

We  use  this  case  to  demonstrate  the  membrane  model  on  a  surface  diffusion-reaction  problem  on 
the  geometry  typical  to  the  Bacillus  type  of  bacteria.  As  discussed  in  Section  4.4  above  the 
JCoBi  software  can  generate  parametric  geometric  shape  of  a  Bacillus  type,  cigar-like  shape,  of  a 
bacterial  membrane.  The  unstructured  surface  triangular  mesh  is  also  automatically  generated  in 
JCoBi  by  simple  geometrical  scaling  based  on  user  specified  geometrical  parameters  of  the 
bacterium  (length  and  radius). 

This  test  case  is  setup  as  a  3D  unsteady  diffusion-reaction  problem  for  two  species,  A  and  B, 
solved  on  the  bacterial  surface  shown  in  Figure  5.4.  The  boundary  conditions  and  the  reaction 
mechanism  are  similar  to  the  test  case  presented  in  Section  5.3  above.  A  simple  one  step 
chemical  reaction  as  follows: 


A-\-  B  '2.B  "1“  B 

Note  that  the  Specie  B  is  both  a  catalyst  for  specie  A  and  a  final  product  of  the  reaction. 

On  the  left  side  of  the  bacterium  the  concentration  of  the  species  B  is  set  to  one  while  the  species 
A  is  set  to  one  on  the  right  side.  The  concentrations  for  both  species  A  and  B  on  the  whole 
membrane  are  initialized  to  zero.  The  geometry,  grid  and  the  concentration  of  specie  B  on  the 
bacterial  surface  at  time  .9  sec.  is  shown  in  Figure  5.4.  The  contours  of  species  B  on  the 
bacterium  membrane  as  a  function  of  time  are  plotted  in  Figure  5.5.  As  before  the  computational 
results  match  our  expected  trends. 
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Figure  5.4:  Geometry  of  a  typical  cigar-like  shape  of  a  bacterial  membrane  and  the 
concentration  profile  of  specie  B  at  t=0.9  sec  . 


t  =  0. 1  sec  t  =  0.7  sec  t  =  0.8  sec  t  =  0.9  sec 


Figure  5.5:  Transient  distributions  of  specie  B  concentration. 

5.5  EGFR  Signal  Transduction  Pathway 

The  epidermal  growth  factor  receptor  (EGFR)  has  been  found  to  be  altered  in  a  variety  of  human 
cancers.  A  number  of  agents  targeting  these  receptors  are  currently  in  clinical  trials  or  are 
already  approved  for  clinical  treatment  [Personal  communication  with  Dr  Lance  Liotta  and  Dr 
Gabriela  Araujo  of  FDA-NCI  Clinical  Proteomics  Program,  Laboratory  of  Pathology,  Center  for 
Cancer  Research,  NCI/NIH,  also  Araujo  et  all  2005].  It  is  generally  accepted  that  better 
understanding  of  EGFR  signaling  pathways  will  be  very  helpful  in  cancer  diagnostics,  designing 
genetic  screening,  co-targeting  additional  proteins  upstream/downstream  of  EGFR  to  improve 
cancer  therapy,  and  to  optimize  caner  treatment. 

EGFR  signaling  pathways  have  been  recently  computationally  analyzed  by  several  investigators 
[Araujo  2005].  We  use  this  test  case  to  validate  CoBi  against  benchmark  results  of  Araujo  2005 
and  to  compare  our  results  to  similar  calculations  performed  with  J-Designer  software  available 
at  CDRC.  J-Designer  is  a  biochemistry  simulation  software  developed  under  DARPA 
BioCOMP  Program  http : //sbw.kgi . edu/ software/ i desi gner.htm 

CoBi  has  been  used  to  reproduce  the  modeling  results  for  the  time  dependent  concentrations  of 
several  species  (RP,  RGS,  RShGS,  and  RPLP)  involved  in  the  EGFR  signal  transduction  cascade 
reported  by  Araujo  et  al.  2005  (Figure  5.6).  Unlike  other  tools,  CoBi  can  be  used  for  detailed 
analysis  of  EGFR  signaling  pathway  combined  with  transient  flow/diffusion  models  and  spatial 
models  of  signal  transduction  in  tissues. 
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Figure  5.6:  Diagram  of  EGFR  signal  transduction  pathway. 

The  EGFR  signaling  pathway  simulations,  reported  by  Araujo  et  al.  2005,  were  repeated  at 
CFDRC  using  both  JDesigner  and  CoBi  codes  and  results  of  time  variations  of  several  species 
concentrations  including  RP,  RGS,  RShGS,  and  RPLP  were  compared.  Both  JDesigner  and 
CoBi  simulations  reproduced  the  results  reported  by  Araujo  [2005]  very  well.  Similar 
computational  times  for  CoBi  and  JDesigner  have  been  observed.  The  simulation  results  of  all 
three  methods  for  species  PRLP  are  shown  in  Figure  5.7. 
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Figure  5. 7;  Time  course  for  specie  PRLP  during  the  EGFR  signaling  pathway  simulations  - 
A.  Araujo  et  al. [2005] ;B.  Using  J Designer;  C.  Using  CoBi. 

5.6  Multicellular  Perfused  Tissue  EGF  Signaling 

In  this  case  we  demonstrate  3D  modeling  of  species  transport  in  the  multicellular  tissue  with 
specific  models  of  convection,  diffusion,  and  reaction  in  perfused  inter-cellular  space 
(extracellular  matrix),  membrane,  and  intracellular  space  (cytosol).  3D  time  dependent  partial 
differential  equations  (PDEs)  were  solved  for  all  species  within  the  intracellular  space  of  the 
cytosol,  in  extracellular  matrix  (ECM)  and  on  the  cell  membrane  surface.  The  species 
conservation  is  expressed  in  the  form: 
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^^  +  u-VC- DV^C  =  R 
dt 


(5.1) 


Where  C  is  the  drug  concentration,  u  is  the  convective  velocity,  D  is  the  drug  diffusivity  (can  be 
anisotropic  based  on  the  cellular  morphology  and  different  in  ECM,  cytosol  and  on  the 
membrane),  R  is  the  generalized  reaction  rate  for  specie  C  as  shown  in  the  reaction  rate 
expressions  table  below.  The  reaction  kinetics  is  similar  to  Lauffenburger  EGFR  model  of 
endocytosis,  Lauffenburher  1999,  Maly  2004,  2004,  and  Wiley  2003.  The  cell-level  model 
simulates  the  time-progression  of  following  species: 

•  Extracellular  ligand  concentration  L, 

•  Surface  receptors  Rs  and  ligand-receptor  complexes  Cs, 

•  Intracellular  ligand  Li,  receptors  Ri,  and  their  complexes  Ci,  and 

•  Degraded  ligands  Ld 

We  use  this  kinetic  model  to  simulate  ligand  depletion  and  endosomal  sorting  after  an  initial 
bolus  of  ligand  into  the  extracellular  medium.  Figure  5.8  presents  the  schematics  of  the  cellular 
endocytosis  and  receptor  recovery  process  and  the  signaling  pathway  used  for  this  model,  Sarkar 
2003.  Figure  5.9  presents  the  kinetic  and  rate  equations  for  EGF  signaling  model  used  in  our 
simulations.  In  most  of  the  previous  simulations  the  process  was  modeled  as  a  well  stirred 
reactor  model.  In  this  project  we  simulate  the  problem  in  a  3D  model  of  the  perfused  tissue. 


43 


Kinetic  Equations 

Rate  Equations 

Extracellular  Ligand  L  df  =  -Vj  -i-  Vj 

v,=kf-R^-k_f:^ 

Surface  Receptor  Rs 

V2  =  k^R^ 

=  -Vi  -  V2  +  Vg 

V3  =  k,C^ 

Surface  Complex  Cs  = 

< 

II 

1 

-0 

II 

Intracellular  Ligand  Li 

3  3  1 

=  -V4  -  V5 

,  so 

II 

Intracellular  Receptor  Ri 

V7  =  k^L. 

d,Ri  =  v^-v,-v, 

v,=V 

Intracellular  Complex  Ci 

=  V3  +  V4  -  V7 

Figure  5.9:  Kinetic  and  rate  equations  for  EGF  signaling  model. 

The  3D  tissue  geometry  was  generated  based  on  mieroscopy  image  of  the  eellular  strueture 
shown  in  Figure  5.10a  below.  To  generate  the  computational  mesh  we  have  used  CFD- 
Micromesh  software  developed  at  CFDRC.  CFD-Micromesh  imports  the  gif  image  of  he 
cellular  tissue,  and  identifies  he  cell  volumes.  The  remaining  part  is  treated  as  the  interstitial 
space.  A  membrane  model  is  specified  at  the  outer  boundaries  of  all  cells.  Finally  an 
unstructured  3D  triangle-based  prismatic  mesh  is  generated,  Figure  5.10b.  Figure  5.10b  also 
presents  the  color  map  of  drug  ligand  concentration  in  the  extracellular  space.  Note  that  the 
ligand  does  not  enter  the  intracellular  space. 


Figure  5.10:  Microscopy  image  of  the  tissue,  3D  unstructured  mesh  with  CoBi  predicted  ligand 
concentration  in  the  intercellular  space  of  the  tissue. 


Figure  5.11  presents  transient  simulation  results  for  two  time  instances  tl  and  t2  >  tl  for  the  drug 
ligand  absorption  in  the  perfused  tissue.  We  present  results  for  three  species:  extracellular 
ligand,  membrane  bound  complex,  and  intracellular  complex.  Note  that  free  ligand  exists  only  in 
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the  extracellular  matrix  (Fig.  5.11  top),  membrane  bound  complexes  Cs  appear  only  on  the 
membrane  (Fig.  5.11  middle),  and  the  intracellular  complex  Ci,  is  present  only  in  the  cytosol 
(Fig  5.11  bottom).  As  observed  in  typical  laboratory  experiments  the  intracellular  complexes 
appear  after  certain  time  delay  we  also  predict  the  membrane  and  intracellular  complexes  with 
time  delay.  Note  that  Ci  is  absent  at  tl  but  is  present  inside  the  cell  cytosol  in  the  proximal  cells 
at  t2  >  tl. 


Ligand 


Complex  Ci 


Complex  Ci 


Ligand 


Complex  Cs 


Time  tl 


Time  t2  >  tl 


Figure  5.11:  Computational  results  for  modeling  EGF  signaling  and  ligand  drug  delivery  to  a 

3D  tissue  structure  with  CoBi. 

5.7  Cellular  Calcium  Oscillations 

A  published,  single  cell  calcium  wave  model  was  used  as  a  test  case  for  CoBi  [McKenzie  1998]. 
This  model  was  selected  because  calcium  waves  have  been  observed  in  single  cells  and  it  is 
believed  that  they  may  play  important  roles  in  signal  transduction  in  several  human  cell  types  in 
both  physiological  and  pathological  conditions.  The  results  reported  for  the  McKenzie  model 
were  reproduced  using  the  CoBi  solver. 

Figure  5.12  shows  a  schematic  of  the  McKenzie  calcium  wave  model.  The  model  uses  a  non- 
physiological  2-D  square  cell  that  is  1  mm  on  a  side.  No  species  are  allowed  to  move  into  or  out 
of  the  cell,  which  comprises  a  mesh  of  40,000  grid  cells.  Calcium  is  normally  stored  in  the 
endoplasmic  reticulum  (ER)  and  is  released  into  the  cytosol  when  IP3  (inositol  triphosphate) 
binds  to  calcium  channels  in  the  ER  membrane.  IP3  is  generated  at  the  cytosolic  side  of  the 
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plasma  membrane  when  certain  types  of  cell  surface  receptors  are  activated.  This  model 
simulates  the  activation  of  receptors  by  using  initial  conditions  in  which  3  separate  areas  of  the 
membrane  have  high  IP3  concentrations  (A,  B,  and  C).  This  model  also  assumes  that  the  ER  is 
uniformly  distributed  within  the  cell  and  ignores  intracellular  compartmentalization. 


Figure  5.12:  Schematic  of  cellular  calcium  wave  model. 

Figure  5.13  shows  the  results  for  various  time  points  during  a  calcium  wave  simulation.  These 
results  are  in  qualitative  agreement  with  the  results  published  for  McKenzie  1998  calcium  wave 
model.  Quantitative  agreement  was  not  possible  because  the  publication  did  not  specify  the  sizes 
and  locations  of  A,  B,  and  C  in  Figure  5.13  and  the  published  figure  was  not  drawn  to  scale. 


Figure  5.13:  Cellular  calcium  wave  simulation. 
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5.8  Tissue  Calcium  Waves 


A  model  of  calcium  wave  propagation  between  adjacent  cells  in  a  tissue  was  reproduced  in  CoBi 
[Sneyd  1998].  This  model  was  selected  because  it  is  believed  that  calcium  waves  may  play  an 
important  role  in  cell-to-cell  communication  within  several  organs  in  humans.  Figure  5.14A  and 
Figure  5.14B  show  schematics  of  the  cell-to-cell  calcium  wave  pathway  and  geometric  mesh  for 
the  model.  The  model  mesh  includes  49  squares,  each  representing  a  biological  cell,  including 
the  plasma  membrane.  Each  cell  square  contains  a  mesh  of  16*16  smaller  squares  that  is  used  to 
model  reactions  and  species  transport  within  each  cell.  Figure  5.14C  shows  a  representative 
result  of  the  CoBi  simulation. 


Figure  5.14:  Tissue  calcium  wave  model  -  A.  Schematic  of  pathway  (Sneyd  1998)  B.  Model 
schematic  (Sneyd  1998)  and  C.  Calcium  concentrations  during  CoBi  simulation. 

The  results  reported  for  the  Sneyd  model  were  reproduced  with  CoBi  except  that  the  fluctuations 
in  cytosolic  calcium  ion  concentration  decayed  more  rapidly  in  our  simulation  and  did  not  decay 
to  baseline  levels  (Figure  5.15).  This  result  appears  to  derive  from  differences  in  initial 
conditions  and  rate  constants  for  calcium  ATPases  that  return  calcium  into  the  ER,  which  were 
not  completely  disclosed  in  the  Sneyd  1998  published  model.  Our  results  do  show  the  damped 
oscillatory  patterns  of  calcium  ion  concentrations  reported  by  Sneyd  1998. 


Figure  5.15:  Comparison  of published  and  CoBi  simulation  results  for  calcium  wave 
transmission  between  cells.  A.  Ca^^  concentrations  (solid  line)  from  Sneyd  and  B.  CoBi. 
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5.9  Morphogenesis  of  Yeast  Cell 


Actin  polymerization  networks  in  yeast  cells  were  simulated  in  CoBi  using  a  model  of  fission 
yeast  cell  division  developed  by  Dr.  John  Tyson  at  Virginia  Polytechnic  University.  The  CoBi 
solver  was  used  to  reproduce  the  time  dependent  changes  in  cell  morphology  for  wild  type  and 
mutant  strains  of  yeast  during  cell  growth  and  division.  Figure  5.16  shows  a  schematic  of  actin 
network  formation  and  breakdown  in  yeast  cells. 


Figure  5.16:  Actin  filament  network  dynamics  [Csikasz-Nagy  2005], 

Protein  U  binds  to  microtubules  and  is  transported  to  the  ends  of  the  cell  where  it  is  released  to 
promote  the  autocatalytic  feedback  loop  for  actin  polymerization.  The  delivery  of  U  to  cell  ends 
by  microtubules  plays  an  important  role  in  localizing  actin  polymerization  to  the  ends  of  the  cell. 
If  the  microtubules  are  disrupted,  the  convection  mechanism  is  stopped  and  U  distributes 
homogeneously  throughout  the  cell.  Figure  5.17  shows  the  results  of  a  simulation  in  CoBi  that 
reproduces  the  bending  of  a  growing  yeast  cell  having  a  mutation  that  causes  a  bent  morphology. 
The  input  data  for  this  test  case  is  presented  in  Appendix  B  to  this  report. 


Figure  5.17:  CoBi  simulation  results  for  yeast  mutants  showing  F-actin  concentrations. 
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5.10  Cellular  Oxygen  and  Energy  Metabolism 
Oxygen  Metabolism  Model 

A  cell  hypoxia/ischemia  simulation  has  been  tested  using  CoBi  for  its  potential  applications  in 
trauma  injury  modeling  for  military  medicine.  A  cell  metabolic  model  based  on  the  cardiac 
myocyte  model  of  Salem  et.  al  [Salem,  2002]  that  includes  responses  to  hypoxic  and  ischemic 
conditions  was  implemented  in  CoBi.  It  models  the  transport  of  glucose,  lactate,  fatty  acids, 
oxygen  (O2)  and  carbon  dioxide  (CO2)  between  the  blood  and  the  cytosolic  compartments.  ATP 
is  produced  from  multiple,  interdependent  pathways  including  glycolysis,  lactate  oxidation,  the 
TCA  cycle,  and  oxidative  phosphorylation. 

The  cellular  oxygen  metabolism  model  comprises  two  compartments,  blood  and  tissue,  with 
passive  gas  diffusion  across  the  compartment  boundaries.  The  diffusive  flux,  J,  of  oxygen, 
across  the  blood-tissue  barrier  is  given  by: 


■  ^b^t,02^t,02) 


(5.2) 


where  is  a  diffusivity  constant  and  is  the  partition  coefficient.  In  general  for  O2  and 
CO2  we  set  =  I .  For  other  species  such  as  free  fatty  acids  and  other  substrates  ^  I . 

The  ordinary  differential  equations  evolving  O2  concentrations  over  time  in  both  blood,  Q  ^2  5 
and  tissue,  Q  ^2  ’  compartments  are  then: 


=  (5.2) 

at 

=  (5.3) 

Equations  for  CO2  have  a  similar  form.  Q  is  the  blood  flow  rate;  the  arterial  or  input 

concentration;  and  F,  are  the  relative  compartmental  volumes  in  blood  and  tissue, 

respectively.  In  this  form,  relative  volumes  are  tuning  parameters  to  calibrate  the  sensitivity  of 
each  compartment  to  perturbations  around  the  steady  state.  The  oxygen  consumption  rate, 
MRq2  ,  is  given  by: 


MR„,  =  a(c,„,  /(*,  +  C,„, ))  RMR  (5.4) 

Here,  RMR  is  the  relative  metabolic  rate,  so  at  the  basal  metabolic  rate,  RMR  =  1 .  The 
parameters  a  and  are  tuning  parameters  and  contain  no  biological  significance.  The  aim 

was  to  introduce  an  upper  bound,  a  ;  to  the  turnover  rate  of  oxygen  because  cellular  oxygen 
turnovers  are  saturable.  After  the  complete  oxidation  of  one  mol  of  glucose,  about  6  mol  of  O2 
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and  CO2  are  consumed  and  produced,  respectively.  Therefore,  it  was  decided  to  set 
MRq2  ~  here  for  simplicity.  In  the  real  cell  glucose  will  be  accompanied  by  fatty  acids, 

lactate,  phosphocreatine  and  ketone  bodies  as  sources  of  energy.  Hence  we  expect  this  simple 
model  to  capture  the  rough  qualitative  features  O2  and  CO2  fluctuations  under  varying  work 
output  (RMR)  levels. 

To  calibrate  the  model  all  equations  were  entered  into  optimization  software.  The  tuning 
parameters  a  ,  k^,  were  solved  for  by  setting  JQ/Jt  =  0  and  dCjdt  =  0 

for  both  species  at  rest  {RMR  =1).  The  four  concentration  initial  conditions,  Q,  and  were 
set  by  organ-specific  experimental  data  [Zhou,  2005]. 

The  reduced  model  is  able  to  capture  qualitative  trends  of  a  normal  response  (Figure  5.18).  In 
response  to  a  doubling  of  the  relative  metabolic  rate  (RMR),  a  new  steady  state  is  reached  in 
which  CO2  output  is  increased  and  blood  and  tissue  oxygenation  is  decreased  slightly. 


Figure  5.18:  Simple  metabolic  model  response  to  a  doubling  of  RMR  from  the  rest  state. 

Energy  Metabolism  Model 

The  present  model  built  in  CoBi  is  an  adaptation  and  improvement  over  a  published  model  for 
cardiac  myocyte  metabolism  [Salem,  2002].  Cardiac  myocytes  are  muscle  cells  and  constantly 
require  large  amounts  of  energy  to  contract  and  pump  blood.  Our  model  consists  of  three 
separate  compartments:  blood,  cytoplasm  and  mitochondria.  This  separation  into  functionally 
separate  compartments  adds  flexibility  and  facilitates  integration  into  the  full  human  model. 
Communication  between  the  three  compartments  occurs  via  passive  or  active  transport 
depending  on  the  species.  The  blood  carries,  oxygen,  lactate,  fatty  acids  and  glucose  to  the  cell 
compartment  for  ATP  synthesis  to  fuel  the  cell  and  carries  carbon  dioxide  and  other  byproducts 
away  from  the  cell.  The  energy  liberated  from  the  conversion  of  ATP  into  ADP  and  AMP  serves 
a  multitude  of  vital  functions  in  the  cell  including  active  transport  across  membranes  and  protein 
synthesis.  Here,  while  some  ATP  will  be  used  in  other  parts  of  the  metabolic  web  (like 
glycogenesis)  the  vast  majority  will  be  converted  into  mechanical  work  for  pumping  blood 
against  a  blood  pressure  differential. 

The  metabolic  scheme  describes  18  key  substrates  located  strategically  in  the  major  nodes  (end 
points  and  branching  points)  of  the  full,  known  metabolic  network  (Figure  5.19). 
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Figure  5.19:  Metabolic  network  diagram  of  the  model  developed  at  CFDRC  (based  on  model  of 

Salem  et.  al.  2002). 


The  model  represents  the  aerobie  and  anaerobic  catabolism  of  glycogen,  fatty  acids  and  lactate 
coming  from  the  blood.  It  also  represents  storage  of  energy  reserves  as  phosphocreatine, 
glycogen  and  triglycerides.  Interconversion  rates  between  substrates  are  quantified  by 
Michaelis-Menten-type  (MM)  kinetics  without  explicit  consideration  of  enzyme  concentration, 
activation,  inhibition  or  other  effects  that  modulate  enzymatic  activity.  Instead,  the  MM 
mathematical  form  is  used  with  the  ratios  of  key  regulatory  factors  in  place  of  substrate 
concentration.  The  regulatory  factors  present  in  the  model  are:  PS  =  [ADP]/[ATP];  RS  = 
[NADH]/[NAD];  CS  =  [CR]/[PC];  and  AF  =  [AC]/[FC].  Our  concentrations  always  refer  to  the 
average  concentration  except  when  explicitly  referred  to  as  blood  compartment  concentrations. 
These  regulatory  factor  ratios  are  known  to  exert  major  regulatory  functions  in  the  cell.  For 
instance,  the  presence  of  high  levels  of  [ATP]  inhibits  the  turnover  rate  of  oxidative 
phosphorylation  while  the  presence  of  [ADP]  stimulates  it.  This  and  other  similar  effects  are 
captured  in  an  approximate,  but  robust  manner.  For  instance,  the  rate  of  conversion  of 
phosphocreatine  to  creatine  is  represented  by 


^PC^CR  ^MAX.PC^CR 


PS 


PS„  +  PS 


c 


PC 


(5.5) 


where  (jfpc^cR  the  flux  or  turnover  rate  (pmoFg-min,  as  all  units  are  on  a  per  gram  of  wet 
tissue  basis),  2.f^Ax,pc^cR  the  maximum  flux  and  PSo  is  the  MM  constant  which  we  set  equal  to 
PS  at  rest.  Cpc  is  the  concentration  of  phosphocreatine  in  the  cytoplasm  (the  ‘c’  subscript  was 
dropped  here  for  simplicity).  All  metabolites  and  their  abbreviations  are  found  in  Table  5.1. 
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Table  5.1:  Metabolite  species  and  their  abbreviations. 


Species 

Abbreviation 

Compartment 

Transport 

Glucose 

GL 

Blood  and  Cytoplasm 

Active  (b  to  c) 

Lactate 

LA 

Blood  and  Cytoplasm 

Active  both  ways 

Fatty  Acids 

FA 

Blood  and  Cytoplasm 

Passive 

Oxygen 

02 

Blood  and  Cytoplasm 

Passive 

Carbon  Dioxide 

C02 

Blood  and  Cytoplasm 

Passive 

Glucose-6-Phosphate 

GP 

Cytoplasm 

None 

Glycogen 

GY 

Cytoplasm 

None 

Triglyceride 

TG 

Cytoplasm 

None 

Pyruvate 

PY 

Cytoplasm 

None 

Acetyl-CoA 

AC 

Cytop.  and  Mitoch. 

None 

Free  CoA 

FC 

Cytop.  and  Mitoch. 

None 

CoA  in  “pool” 

CoApool 

Cytoplasm 

None 

Nicotinamide  Adenine  Dinucleotide 
(oxidized) 

NAD 

Cytop.  and  Mitoch. 

None 

Nicotinamide  Adenine  Dinucleotide 
(reduced) 

NADH 

Cytop.  and  Mitoch. 

None 

Adenosine  Diphosphate 

ADP 

Cytop.  and  Mitoch. 

None 

Adenosine  Triphosphate 

ATP 

Cytop.  and  Mitoch. 

None 

Phosphocreatine 

PC 

Cytoplasm 

None 

Creatine 

CR 

Cytoplasm 

None 

When  more  than  one  regulatory  factor  influences  the  rate  a  “ping-pong”-type  kinetic  formula  is 
used.  For  instance,  the  conversion  of  pyruvate  to  acetyl-CoA  is 


^PY^AC  ^MAX.PY^AC 


HRS. 


1  +  ^^ — -  + 

l/AF  l/RS 


(5.6) 


The  use  of  these  forms  correctly  reduce  reaction  rates  to  zero  as  the  concentrations  of  reactants 
reach  zero.  This  is  an  improvement  over  previous  models,  which  suffered  from  non- 
physiological  results  such  as  negative  concentrations  under  all  but  a  very  narrow  set  of 
physiological  input  parameters  for  which  the  model  was  designed  to  mimic.  The  application  of 
these  robust  and  simple  rates  allows  this  model  to  be  used  in  a  much  wider  range  of 
physiological  states  without  retuning  the  rate  constants.  The  full  compliment  of  reaction  rate 
equations  can  be  found  in  Appendix  B. 


All  metabolite  concentrations  are  evolved  in  time  by  integration  of  ordinary  differential 
equations  (ODEs)  expressing  mass  conservation.  All  ODEs  have  the  general  form: 

volume  X  —  (concentration)  =  production  -  utilization  +  flux  in  -  flux  out. 
dt 


In  the  cytoplasm,  this  is  written  as: 

n 


dC, 

dt 


a, 


k=\ 


z 

k=\ 


a 


jk 


tF  or  P 


jtort-  _  jt 

^  b^cj  c->-bJ 


(5.7) 


and  similarly  to  the  mitochondrial  compartment,  m.  The  above  equation  is  expressed  in  the 
symbols  that  will  be  used  in  the  remainder  of  this  text.  Here,  the  volume,  Vj ,  is  the  effective 


compartmental  volume  and  a  tuning  parameter  for  adjusting  the  sensitivity  of  the  metabolite  to 
perturbations  to  the  independent  variables.  In  essence,  it  adjusts  the  timescale  (l/Vj )  of  the 
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model.  The  production  and  utilization  terms,  the  two  summations,  respectively,  are  just  the 
algebraic  sum  of  the  rates  of  conversion  of  species  ‘j’  to  species  ‘k’,  or  ,  multiplied  by  their 

corresponding  stoichiometric  coefficients,  a^.j  (all  stoichiometric  equations  are  given  below  in 
Table  5.2).  The  fluxes  in  and  out  of  the  cytoplasm  of  species  ‘j’  are  the  mass  transport  rates 
JbXlj  '^c-Zyj-’  respectively.  These  fluxes  may  be  passive  (superscript-P),  or  active 

(superscript-F),  depending  on  the  species.  In  the  blood  compartment  there  are  no  production  or 
utilization  terms.  The  general  form  is  written: 

Vb^  =  Q(C^j-C,^)-j[:Z  (5-8) 

here,  Q  represents  the  coronary  blood  flow  rate,  which  is  set  to  1.0  mL/min  at  rest.  The 
concentrations  of  arterial  and  mixed  venous  blood  are  C^j  and  C^j,  respectively.  When 

calculating  the  transport  rates,  the  relevant  quantities  are  not  the  arterial  or  venous  blood 
concentrations,  but  a  mixture  of  the  two.  To  remedy  this,  an  average  blood  concentration,  is 
written  which  assumes  a  linear  combination  of  arterial  and  venous  bloods  at  a  proportion  set  by 
the  mixing  fraction,  F=  0.33: 

C,,=FQ,+(1-F)C,,.  (5.9) 


Table  5.2:  Summarized  stoichiometric  reactions. 


Reaction _ Equation 

CYTOSOL 


Glycolysis  1 
Glycolysis  2 
Glycogenesis 
Glycogenolysis 
Pyruvate  Oxidation 
Pyruvate  Reduction 
ATP  hydrolysis 
Palmitate  Oxidation 
Triglyceride  Production 
Triglyceride  Oxidation 

Phosphocreatine  Synthesis 

Phosphocreatine 

Breakdown 


Glucose  +  ATP^  Glucose-6-phosphate  +  ADP  + 

Glucose-6-phosphate  +  2  Pi  +  3  ADP  +  2  NAD^^2  Pyruvate  +  3  ATP  +  2  NADH  +  +  2  H2O 

Glucose-6-phosphate  +  ATP  +  Glycogen  (n)  +  H2O  ^  Glycogen  (n+1)  +  ADP  +  2  Pi 
Glycogen  (n+1)  +  Pi  ^  Glycogen  (n)  +  Glucose-6-phosphate 
Pyruvate  +  CoA  +  NAD^  ^  Acetyl-CoA  +  CO2  +  NADH 
Pyruvate  +  NADH  +  H^  ^  Lactate  +  NAD^ 

ATP  ^  ADP  +  Pi  +  Energy 

Palmitate  +  8  CoA  +  35/3  NAD^  +  2  ATP  +  7  H2O  8  acetyl-CoA  +  35/3  NADH  +  7  H^  +  2  ADP 
+  Pi 

3  Palmitate  +  Glycerol  +  2  ATP  ^  Triglyceride  +  2  ADP  +  2  Pi 
Triglyceride  ^  3  Palmitate  +  Glycerol 
Creatine  +  ATP  ^  Phosphocreatine  +  ADP  +  H^ 

Phosphocreatine  +  ADP  +  H^  ^  Creatine  +  ATP 


MITOCHONDRIA 


Oxidative  Phosphorylation  NADH  +  H^  +  3  ADP  +  3  Pi  +  O2  ^  NAD^+  4  H2O  +  3  ATP 

Citric  Acid  Cycle  Acetyl-CoA  +  1 1/3  NAD^  +  ADP  +  Pi  +  2  H2O  1 1/3  NADH  +  2  H^  +  CoA  +  ATP  +  2  CO2 
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As  mentioned  previously,  the  mass  transport  rates  can  be  either  passive  or  active.  The  passive 
rate  is  set  by 

(5.10) 

where  and  y  are  adjustable  parameters  controlling  the  equilibrium  set  point  and  gain 

of  the  flux.  Notice  also  that  similar  to  a  difflisivity  constant.  Blood  gasses  and  fatty 

acids  are  transported  passively.  In  the  case  of  active  transport,  the  active  sites  available  to  carry 
molecules  across  membrane  barriers  are  finite  and  therefore  saturable.  Hence,  we  employ  an 
approximate  rate  of  saturable  form,  which  resembles  a  MM  rate: 


b^cj 


^ b^cj  ^b,j 


(5.11) 


Like  in  the  case  of  the  reaction  rate  equations,  we  set  the  constant  in  the  denominator,  j . 
equal  to  the  value  for  C^j  at  rest. 


Coupling  to  Work  Rate 

Oxidations  of  fuels  is  controlled  such  that  the  ATP  requirements  for  a  variety  of  vital  process 
within  the  cell  are  met.  Therefore,  even  though  energy  conversion  processes  within  the  cell  are 
typically  efficient  (7«0.5),  not  all  of  the  free  energy  not  converted  into  heat  goes  into 
mechanical  work  production.  Here  we  assume  that  only  a  small  fraction  =  Q.\  contributes  to 
the  mechanical  pumping  of  blood  throughout  the  body.  Also,  in  calculating  how  much  energy  is 
available  from  the  hydrolysis  of  ATP  to  ADP  we  have  neglected  AMP  forms  and  written 

^^ATP^ADP  ~  ATP ^ ADP  +-^Tln(Cp.  •  PS^  (5-12) 

This  expression  also  neglects  other  factors,  such  as  pH  and  [Mg^^],  which  also  influence  the 
available  free  energy.  Furthermore,  because  our  model  does  not  treat  free  phosphate  (Pi) 
fluctuations,  we  fix  Cp,  to  the  concentration  at  rest,  letting  only  PS  vary  over  time.  The  standard 

free  energy  of  ATP  hydrolysis  is  set  to  =-31.5  kJ/mol. 


The  rate  of  ATP  hydrolysis  is 


ATP^ADP 


W _ l/PS 

—  pAG ^pp^^jjp  •  M  K^pp^yiDP  yPS 


(5.13) 


where  W  is  the  work  rate  (power)  demand  for  pumping  blood,  which  is  equal  to  the  product 
QAP ,  or  the  cardiac  output  times  the  pressure  load  across  the  heart.  M  is  the  mass  of  wet  heart 
tissue  since  all  expressions  are  written  on  this  basis.  The  second  factor  in  the  equation  was  added 
ad  hoc  to  force  the  conversion  rate  to  vanish  as  ATP  is  depleted.  In  this  sense,  the  first  term 
represents  the  maximum  work  rate  demands  by  the  heart  (similar  to  Fmax),  while  the  second 
term  modulates  the  actual  energy  delivery  according  to  the  availability  of  ATP.  Setting 
Katp^adp  «  o’ conditions  near  the  rest  state,  allows  the  model  to  supply  almost  the  full 
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amount  of  energy  demand.  As  the  heart  works  harder,  the  supply-to-demand  ratio  will  drop  as 
ATP  becomes  depleted. 

Model  Optimization 

Because  the  general  topology  of  the  metabolic  network  is  conserved  among  many  different  cell 
types,  the  current  model  could  be  calibrated,  or  optimized,  to  represent  different  tissues  with 
perhaps  only  minor  modifications.  Where  other  metabolic  processes  are  important,  these  extra 
pathways  could  be  added  to  the  basic  model  and  recalibrated  accordingly. 

Here  the  rate  constants  were  optimized  at  steady  state  rest  conditions  to  mimic  the  heart  muscle 
(myocytes).  The  strategy  was  to  obtain  a  constant  concentration  of  all  metabolites  for  rest  levels 
of  arterial,  blood,  cytoplasm  and  mitochondria  concentrations  by  tuning  the  reaction  and 
transport  rate  constants.  The  steady  state  constraint  is  simply: 

dC, 

— ^  =  0  (5.14) 

dt 

for  all  concentrations  in  the  blood,  cytoplasm  and  mitochondria  compartments.  The  initial 
conditions  at  rest  are  also  known  (or  estimated)  and  fixed.  In  addition,  the  MM  constants  {Km's) 
were  set  to  their  rest  values  as  mentioned  previously.  With  this,  the  only  variables  in  the 
optimization  were:  all  ,  M ,  ,  and  the  transport  coefficients  and  j . 

Finally,  in  the  heart,  about  60  to  90%  of  the  acetyl-CoA  comes  from  the  breakdown  of  fatty 
acids.  This  was  added  as  an  extra  constraint  during  the  fine-tuning  of  model  parameters: 

- tf±iA£ - =  0.67  at  rest.  (5.15) 

^FA-^AC  ^PY^AC 

Optimization  was  accomplished  by  a  Newton  search  algorithm  with  an  initial  guess  estimated 
from  There  were  26  variables  and  24  equations  in  the  final  optimization.  The  two  degrees  of 
freedom  in  the  model  were  used  to  perform  a  limited  amount  of  manual  fine-tuning  on  the 
sensitivity  of  the  response  to  better  match  available  experimental  data  (see  Model  Validations). 
The  timescale  of  the  responses  (the  time  to  reach  steady  state  following  a  perturbation)  were 
adjusted  by  varying  the  effective  compartmental  volumes  ( Vj )  until  it  approximately  matched 

experimental  data  on  the  evolution  of  several  metabolites.  The  relative  values  of  the  volumes, 
however,  were  fixed  according  to  published  heart  data  *  to  87%  for  heart  and  13%  for  blood. 

Model  Validations 

The  predictive  power  of  our  model  was  tested  against  available  experimental  data  on  dog,  swine 
and  pig  models  [Pantely,  1990,  1991;  Stanley  1992].  In  these  studies,  blood  flow  through  the  left 
anterior  descending  (LAD)  coronary  artery  was  reduced  to  prescribed  fractions  of  the  control 
subject  for  1-hour  in  anesthetized  animals.  These  ischemia  levels  were  set  to  30%  in  [Pantely, 
1990];  50%  in  [Pantely,  1991];  and  60%  in  [Stanley  1992]. 
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Simulations  lasted  65  minutes.  After  5  minutes  of  evolution  at  rest  levels,  a  step  change  was 
applied  to  the  coronary  blood  flow  parameter  Q  to  the  desired  level  of  ischemia.  The  width  of  the 
step  change  was  1  minute  but  alterations  of  width  did  not  affect  the  resulting  fluctuations 
significantly.  In  general,  simulations  yielded  well-behaved  functions  except,  in  a  few  cases,  for 
very  early  abrupt  fluctuations  following  immediately  after  the  change  in  Q.  The  product  of  heart 
rate  and  left  ventricular  systolic  pressure  remained  unchanged  during  several  ischemia 
experiments  of  varying  intensities.  In  general,  while  the  heart  rate  increased  slightly,  the  blood 
pressure  decreased  accordingly  to  keep  the  product  nearly  constant.  From  this  it  is  speculated 
that  the  mechanical  power  output  from  the  heart  does  not  vary  appreciably  during  ischemic 
conditions.  In  our  model  we  reflect  this  by  letting  the  mechanical  power  demand  remain  constant 
throughout  the  entire  simulation. 

All  1-hour  ischemia  simulations  ran  in  a  matter  of  just  a  few  seconds  using  a  stiff  integral 
method  with  relative  and  absolute  tolerances  of  10'^. 

According  to  Stanley  et.  al.  [Stanley  1992],  even  moderate  ischemia  produces  little  to  no  change 
in  glucose  uptake  even  though  it  increases  glucose  extraction.  Uptake  is  measured  as  the  product 
of  flow  and  arterial-to- venous  concentration  difference: 

Also,  FA  uptake  decreases  while  the  heart  switches  from  a  net  consumer  to  a  net  producer  of  LA. 
Our  model  was  able  to  reproduce  all  of  these  uptake  patterns  (Figure  5.20). 


Figure  5.20:  Effect  of  GL,  FA  and  LA  uptake  (bottom  graph)  to  30%  ischemia. 

ATP  and  PC  levels  were  also  monitored  in  these  in  vivo  experiments.  Figure  5.21  shows  both 
data  and  simulation  of  these  fluctuations. 
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Figure  5.21:  ATP  and  PC  evolutions  compared.  Experimental  data  points  are  plotted  with 
simulation  curves  for  30%  and  50%  ischemia. 

Comparisons  between  experiment  and  our  simulations  are  summarized  in  Table  5.3.  Because  the 
control  metabolite  concentrations  differed  slightly  between  experiment  and  simulation,  all  values 
were  normalized  with  the  control  (initial)  concentrations  {Cj  =  Cj  jCj^ ).  ATP  and  PC  data  are 
mean  values  of  subendocardium,  midmyocardium  and  subepicardium  levels. 

Table  5.3:  Comparison  of  experiment  to  simulation  during  1-hour  ischemia. 


Metabolite 

5-10  mins 

30%  ischemia 

5  mins 

50%  ischemia 

Experiment 

Simulation 

Experiment 

Simulation 

ATP 

0.70±0.26 

0.83 

0.65±0.10 

0.67 

Phosphocreatine 

0.68±0.15 

0.64 

0.48±0.21 

0.40 

Lactate 

2.26±1.1 

2.14 

02  consumption 

0.77±0.16 

0.92 

In  these  experiments,  the  authors  have  reported  recovery  of  ATP,  PC  and  other  metabolite  levels 
after  60  minutes  of  continued  mild  to  moderate  ischemia.  Arai  et.  al.  speculates  that  active  down 
regulation  could  have  decreased  the  demand  for  ATP  thus  putting  the  heart  in  a  state  of  “partial 
hibernation”.  This  recovery  could  be  a  result  of  intracellular  pH  improving  on  a  similar  time 
course  to  the  myocardial  lactate  content.  In  this  picture,  the  free  energy  available  per  ATP  would 
have  recovered  towards  preischemic  levels.  The  authors  also  point  out  that  pH  alone  cannot 
account  for  full  recovery  and  further  speculate  the  involvement  of  free  phosphate  (Pi),  altered 
Ca^^  transient  size  or  altered  myofibril-Ca^^  interactions  in  the  recovery  process. 

Yet  another  possible  cause  for  recovery  not  pointed  out  by  the  authors  is  the  influence  of  all 
other  tissues  and  organs  in  the  anesthetized  animal  model  via  venous  return.  Since  a  fully 
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coupled  organism  will  differ  from  an  isolated,  heart-only  model  with  constant  arterial  blood 
metabolite  levels,  a  realistic  comparison  to  experimental  models  has  to  be  limited  to  the  earliest 
effects  of  acute  ischemia.  That  is,  before  the  onset  of  possible  compensatory  effects  by  other 
parts  of  the  organism  including  hormonal  regulation.  Thus,  the  comparison  is  limited  to  only  the 
first  5  to  10  minutes,  when  the  first  measurements  were  taken.  Nevertheless,  because  the  model 
is  coupled  to  the  mechanical  power  output  demands  of  the  heart  we  can  easily  test  Aral’s 
proposition  of  reduced  energy  demand  by  asking  how  large  a  decrease  in  the  power  output  could 
have  led  to  a  recovery  of  metabolites  to  preischemic  levels.  It  was  found  that  a  25%  increase  in 
the  available  free  energy  around  the  middle  of  the  experiment  would  have  been  sufficient  to 
return  metabolite  concentrations  to  preischemic  levels  after  the  onset  of  a  30%  ischemia. 

Figure  5.22  shows  a  few  more  results  from  the  30%  ischemia  simulation.  It  is  known  that  in  the 
heart  about  60  to  90%  of  the  acetyl-CoA  comes  from  the  breakdown  of  fatty  acids  (FA  ^  AC) 
with  the  remainder  coming  mostly  from  pyruvate  breakdown  (PY  ^  AC),  both  before  and  after 
ischemia.  This  was  captured  in  our  simulations.  At  rest,  the  FA  to  AC  pathway  represents  67% 
of  the  total  AC  production  (Figure  5.232,  left  panel).  After  the  onset  of  ischemia  this  decreases 
temporarily  and  then  approaches  the  previous  fraction  after  1  hour  of  continued  ischemia.  This  is 
also  reflected  in  the  fluctuations  of  NAD  and  NADH  (Figure  5.232,  right  panel).  Here,  ischemia 
causes  a  decrease  in  the  redox  state.  This  result  is  consistent  with  a  recent  paper  [Zhou,  2005]  but 
contradicts  an  earlier  paper  [Salem,  2002]  by  the  same  group. 


Figure  5.22:  Effects  of  30%  ischemia  on  AC  production  and  redox  state. 

The  coupling  to  mechanical  work,  a  unique  feature  of  this  model,  deserves  special  attention. 
Figure  5.23  shows  the  effect  of  50%  ischemia  on  the  various  parameters  controlling  mechanical 
energy  consumption. 
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Figure  5.23:  Effect  of  50%  ischemia  on  ATP  hydrolysis  rate,  max  ATP  hydro,  rate  and  available 

free  energy  (inset). 

The  graph  shows  that  during  ischemic  conditions  the  maximum  rate  of  ATP  hydrolysis  ( ) 
actually  increases  due  to  the  smaller  free  energy  available  from  the  reaction.  The  rate  of 
hydrolysis,  on  the  other  hand,  decreases  due  to  inadequate  perfusion. 

In  conclusion,  even  thought  our  model  falls  within  the  (rather  large)  error  bounds  of 
experimental  dog,  pig  and  swine  models,  the  mathematical  model  remains  of  qualitative 
precision  only.  This  is  because  of  the  several  shortcomings  of  the  model  including  vastly 
simplified  metabolic  topology  utilized  in  its  construction  and  the  limited  experimental  data 
available  for  model  validation.  However,  the  model  is  fairly  robust  and  will  give  reasonable 
qualitative  results  under  most  physiological  conditions. 


Hypoxia  Simulation 

A  hypoxia  simulation  was  also  performed  during  model  validation  (Figure  5.24).  Even  though 
these  results  were  not  compared  against  experiment  at  this  point,  they  conform  to  expectation. 
Briefly,  they  resemble  ischemia  of  a  more  severe  intensity.  Like  ischemia,  the  heart  switches  to  a 
net  producer  of  lactate  due  to  the  limiting  oxidative  capacity  of  oxidative  phosphorylation.  Fatty 
acid,  oxygen  and  carbon  dioxide  uptakes  are  all  down-regulated  to  cope  with  the  decreased 
oxidative  potential  of  the  cell. 
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Figure  5.24:  Effects  of  40yo  hypoxia  on  the  uptake  of  several  metabolites  and  blood  gasses. 

Model  Stability 

Our  model  was  able  to  display  a  well-behaved  response  (no  oscillations)  for  step  input  changes 
of  as  fast  as  0.1  minutes.  To  demonstrate  stability  we  have  simulated  an  extreme  near  total 
ischemia  (99.9%)  for  1  hour  (Figure  5.25).  Most  physiological  changes  in  the  human  body  are 
relatively  milder  and  slower  changing,  however.  Thus  it  is  expected  that  the  current  model 
remains  stable  after  incorporation  into  the  larger-scale  human  model. 


Figure  5.25:  99.9%  ischemia  simulation. 
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5.11  Perfusion  of  a  Cell  in  an  Organ  and  Cell  Metabolism 


A  realistic  spatio-temporal  metabolism  model  was  developed  in  CoBi.  The  CoBi  metabolism 
model  consists  of  a  3D  spatially  meshed  model,  using  realistic  single  cell  dimensions  and  an 
unstructured  grid.  The  simulation  was  transient,  and  utilized  CoBi’s  flow,  membrane,  species, 
and  reaction  modules.  The  geometry  consisted  of  capillary  blood  flow  around  a  circular  cell, 
which  included  the  cytosol,  20+  mitochondria,  and  an  unmeshed  nucleus.  Membrane  structures 
were  present  between  the  blood-cytosol  and  the  cytosol-mitochondria  interfaces.  Present  through 
the  geometry  were  18  species.  In  the  blood  flow,  species  were  transported  through  both 
convective  and  passive  transport.  Throughout  the  cytosol  and  mitochondria  species  were 
transported  via  passive  and  active  transport.  Active  transport  was  assumed  to  be  a  saturable  two- 
step  process  where  transporter  sites  in  the  membrane  were  bound  to  the  molecule  before  being 
transported.  This  treatment  also  gave  Michaelis-Menten  kinetics  for  active  transport.  Passive 
membrane  transport  was  governed  through  a  conjugate  gradient  equation.  In  the  mitochondria, 
there  were  7  species  reactions,  while  17  species  reactions  were  occurring  in  the  cytosol.  These 
reactions  in  conjunction  with  the  transport  of  the  species  resulted  in  quasi-steady  state  conditions 
for  the  cell  similar  to  those  seen  in  the  three-compartment  model  (Figure  5.26  and  Figure  5.27). 
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Figure  5.26:  Snapshots  of  the  cell  model.  (Top  Left)  The  unstructured  mesh  grid  used  for 
simulations.  (Bottom  Left)  Oxygen  distribution.  (Top  and  Bottom  Right)  Carbon  dioxide 
distributions.  Capillary  blood  flow  is  shown  circling  both  sides  of  the  cell  and  entering  from  the 
bottom  part  of  the  figure.  Species  are  transported  in  and  out  of  the  cell  through  the  cell 
membrane  and  then  into  the  mitochondria  where  energy  is  converted  into  ATP  and  CO 2. 
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Figure  5.27:  A  series-in-time  showing  the  progression  of  metabolism  in  a  cell.  Top  figures  show 
CO2  concentration  evolution  throughout  the  cell  for  two  different  metabolic  rates  (MR  =  9  and 
90).  The  time  stamps  are  indicated  at  the  top  of  the  figure.  Bottom  figures  show  the  O2 
distribution  for  the  same  simulation  and  time  points.  Notice  the  increased  O2  consumption  and 
CO 2  generation  when  energy  demands  on  the  cell  are  increased  10-fold. 
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6. 


CONCLUSIONS  AND  RECOMMENDATIONS 


6.1  Summary  and  Conclusions 

This  report  presents  the  results  of  a  two-year  project  performed  by  CFD  Research  Corporation 
under  the  DARPA/IPTO  BioCOMP  Program.  The  BioCOMP  Program  objective  was  to  develop 
advanced  simulation  framework  called  Simulation  Program  for  Intra-Cellular  Processes  (Bio- 
SPICE),  for  comprehensive  simulations  of  cell  biology  processes.  The  role  of  the  Bio-SPICE 
framework  is  to  enable  the  construction  of  sophisticated  models  of  intracellular  processes  that 
can  be  used  to  predict  and  control  the  behavior  of  living  cells. 

Most  of  the  cell  biology  models  typically  reported  the  open  literature  are  formulated  as  single 
compartment  (well  stirred)  reaction  models.  Our  approach  in  this  project  was  to  formulate  them 
as  ID,  2D  and  even  3D  transient  (spatiotemporal)  models  and  to  perform  parametric  simulations 
to  assess  model  performance  in  the  multi  dimensional  form.  Spatiotemporal  modeling  of  cell 
biology  is  performed  on  a  numerical  grid.  In  this  project  we  have  explored  two  ways  to  generate 
computational  grids  for  biological  cells:  image  based  mesh  generation,  and  meshing  idealized 
geometries  representative  of  biological  ells  e.g.  sphere  and  rod.  In  this  project  we  have  used 
CFD-Micromesh  software  it  for  image  processing,  geometry  construction,  and  for  generation  of 
2D  and  quasi  3D  (extruded  2D)  meshes  for  simulation  of  cellular  biophysics  and  biochemistry. 
In  the  CoBi  software  environment  the  mesh  is  either  imported  as  a  file  if  the  geometry  is 
complex  or  automatically  generated  in  the  JCoBi  tool  for  simple  geometric  shapes  such  as 
square,  cube,  rod,  sphere,  ellipsoid,  and  bacteria-like  shape. 

A  critical  review  of  published  spatiotemporal  models  describing  biochemistry  and  biophysics  of 
prokaryotic  and  eukaryotic  cells  was  conducted.  The  objective  was  to  select  the  models  with 
important  applications  in  military  medicine.  Several  problems  were  selected  for  parametric 
simulations  including: 

•  Cellular  chemosensing 

•  Cell  cycle  biochemistry  model 

•  Cell  cycle  and  division  model 

•  Polarization  and  oscillations  during  E  Coli  cell  division 

•  Signaling  pathways  during  cacterial  chemotaxis 

•  Multivalent  binding 

•  Neutrophil  chemotaxis 

•  E  Coli  chemotaxis 

The  main  objective  of  the  BioCOMP  Program  was  to  develop  a  self-contained  software 
framework  solving  partial  differential  equations  governing  flow,  diffusion,  and  biochemistry 
applicable  for  spatiotemporal  modeling  of  cell  and  organ  biology  problems.  During  a  two  year 
period  CFDRC  has  developed  CoBi  to  simulate  complex  cell  and  organ  biology  problems.  The 
code  was  written  in  modem  C++  programming  language  using  state  of  the  art  numerical 
techniques  for  solving  large  systems  of  partial  differential  equations  (PDEs)  and  CFDRC 
multiphysics  modeling  expertise.  A  graphical  user  interface  (GUI)  for  CoBi,  JCoBi,  was  written 
in  Java  language  and  interactive  3D  graphics.  CoBi  has  been  designed  to  interact  with  the  other 
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Bio-SPICE  software  tools  using  SBML  as  a  standard  interface.  The  prototype  CoBi  software 
developed  in  this  project  was  delivered  in  source  to  DARPA  (SRI  Bio-SPICE  site)  in  February 
2005  and  the  final  code  along  with  the  code  documentation  was  delivered  to  AFRL/IFTC  at  the 
end  of  the  project. 

During  this  project  CoBi  code  has  been  tested  on  several  fluid  flow,  diffusion,  and  chemistry 
problems  and  results  have  been  compared  to  CFD-ACE-i-  and  to  published  benchmark  solutions. 
The  code  showed  very  good  accuracy  and  numerical  robustness  in  solving  3D  problems.  In  the 
current  project  some  of  the  cell  biology  problems  studied  in  our  previous  BioCOMP  Program 
presented  in  Appendix  A  were  reproduced. 

The  major  achievement  of  this  project  was  the  development  of  a  3D  volumetric  modeling  tool, 
CoBi,  designed  for  spatiotemporal  modeling  of  cell  biology  by  coupling  cell  biophysics  and 
biochemistry  equations.  A  generalized  biochemistry  solver  module  for  modeling  cell  signaling 
and  metabolic  pathways  problems  including  extracellular,  intracellular,  and  membrane  bound 
biochemical  reactions  was  developed.  This  report  describes  the  mathematical  formulation  of 
cellular  biophysics  and  biochemistry  used  in  the  code  development,  numerical  methods  for 
solving  large  scale  systems  of  PDE  equations,  and  CoBi  software  structure.  The  code  has  been 
applied  to  a  number  of  cell  biology  problems.  This  report  describes  details  of  model  description 
and  presents  simulation  results  for  the  flowing  cell  biology  problems: 

•  Bacillus  Subtilis  Min  protein  dynamics 

•  Autocrine  EGFR  signaling 

•  Multi-cellular  perfused  tissue  EOF  signaling 

•  Cellular  calcium  oscillations 

•  Tissue  calcium  waves 

•  Cellular  oxygen  metabolism 

•  Morphogenesis  of  yeast  cell 

•  Perfusion  of  a  cell  in  an  organ  and  cell  metabolism 

At  the  end  of  the  project  CFDRC  has  prepared  CoBi  documentation:  JCoBi  Users  Manual  and 
four  CoBi  tutorials  and  delivered  to  AFRL/IFTC,  Rome  NY,  along  with  the  CoBi  source 
software.  They  are  available  to  other  DoD  Laboratories. 

It  is  felt  that  all  of  the  major  objectives  of  this  project  have  been  successfully  accomplished.  The 
CoBi  code  has  been  developed,  validated,  and  demonstrated  on  several  cell  biology  and  organ 
physiology  problems.  The  code  has  been  used  as  the  starting  point  for  several  DARPA  and  DoD 
projects.  CFDRC  is  using  CoBi  as  the  enabling  technology  for  projects  related  to  military 
medicine,  personnel  protection,  and  bio  defense. 

The  CoBi  code  delivered  at  the  end  of  this  project  will  be  available  to  other  US  government 
organizations  and  to  collaborating  academic  institutions  without  license  fees.  Technical  support 
for  the  effective  use  of  the  CoBi  code  is  available  from  CFDRC. 
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6.2  Future  Developments  and  Recommendations 

This  project  has  established  excellent  foundations  for  future  computational  medicine  and 
biology  projects  (CMB)  at  CFDRC  for  military  and  civilian  medical  applications.  After  two 
years  of  this  project  CoBi  has  matured  to  the  point  that  it  can  be  used  for  solving  problems  not 
only  in  cell  biology  but  also  organ  biology  and  whole  body  (organism)  physiology.  During 
2005-2006  CFDRC  has  invested  significant  IR&D  resources  in  CoBi  code  development, 
integration,  and  testing.  The  code  is  being  used  as  a  starting  point  in  several  projects  with  DoD 
focused  on  biotechnology,  nanotechnology,  and  military  medicine.  As  a  result  several  new 
capabilities  are  being  implemented  in  CoBi  including:  flow  in  porous  media,  conjugate  heat 
transfer,  electrostatics,  electrochemistry,  stress/deformation  and  structures  dynamics,  moving 
body  simulation,  and.  Some  of  the  current  projects  that  utilize  the  CoBi  software  capability  and 
contribute  to  further  development  and  applications  are: 

•  AFRL/IFTC  SBIR  Phase  II  project  on  “Bio-CAD  Tools  for  DNA  Computing”  adapting 
the  CoBi  code  to  simulate  a  DNA  based  bio  computer,  and  a  DNA  based  self  assembly  of 
nano  devices  for  next  generation  of  computers  and  sensors, 

•  DARPA/DSO  BAA  Surviving  Blood  Loss  (SBL)  in  collaboration  with  University  of 
Alabama  Birmingham  and  US  Army  WRAIR  on  “Enhanced  Survival  following  Trauma- 
Hemorrhage  and  Extended  Period  of  Hypotension  Using  Hormonal  and  Metabolic 
Support”.  CFDRC  is  using  CoBi  to  develop  whole-body  model  of  human 
cardiopulmonary  circulation,  respiration,  metabolism  and  hemorrhage.  In  an  ongoing 
project  we  are  also  developing  rat  and  pig  models  to  study  novel  pharmacological 
treatment  and  resuscitation  strategies, 

•  DARPA/DSO  SBIR  Phase  II  on  “Multiscale  Modeling  Tools  for  Blast  Wave  Traumatic 
Brain  Injury”.  In  this  project  CoBi  is  being  adapted  for  modeling  tissue  biomechanics 
(stress/deformation)  during  shock  wave  impact  and  penetration  through  human  lungs  and 
brain, 

•  US  Army  Aeromedical  Research  Laboratory,  Ft  Rucker  AL  Phase  II  project  on 
“Image/Model  Based  System  for  Optimized  Helmet  Design”.  Here  CoBi  is  being  adapted 
to  simulate  kinetic  and  blast  impact  on  a  human  head  with  a  protective  helmet, 

•  AFRL/HEPA  Wright  Laboratories  Dayton  OH  on  “An  Integrated  Modeling  Framework 
for  Predictive  Human  Performance”.  In  this  project  we  will  integrate  CoBi  software  tools 
with  bio-dynamics  and  anthropologic  databases  of  human  (aircraft  pilot)  body  injury, 

•  US  Army  WRAIR  Phase  II  STTR  Project  with  UC  Irvine  on  ’’Portable  Near  Infrared 
Detection  of  Blast  Lung  Injury”  CFDRC  is  developing  integrated  hardware  NIRS  sensor 
and  CoBi  software  for  model  based  detection  of  lung  injury, 

•  NIH  phase  II  SBIR  project  with  Vanderbilt  University  on  “Multiscale  Model  of  Tumor 
Angiogenesis”.  CoBi  is  being  adopted  for  modeling  tumor  perfusion  and  role  of  anti 
angiogenesis  drugs  in  optimized  chemotherapy. 

CFDRC  is  also  involved  in  several  commercial  applications  for  civilian  medicine  and  for  model 
based  development  of  novel  medical  instruments  and  sensors. 
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APPENDIX  A 


Summary  of  CFDRC  Subcontract  Project  with  Lawrence  Berkeley 
Laboratory  under  DARPA  BAA  BioCOMP  Project  on 


An  Experimental  Test  Suite  for  Development  of  the 

BioSPICE  Kernel 


Lawrence  Berkeley  Laboratory  Principal  Investigator: 

Dr.  Adam  Arkin 


CFDRC  Subcontract  Number:  6513845 


Background 

CFD  Research  Corporation  (CFDRC)  has  been  involved  in  the  DARPA/IPTO  BAA  01-26 
BioCOMP  Program  between  January  2002  and  April  2006.  The  objective  of  the  DARPA 
BioCOMP  Program  was  to  develop  a  computational  framework  that  enables  the  construction  of 
sophisticated  models  of  intracellular  processes  that  can  be  used  to  predict  and  control  the 
behavior  of  living  cells.  One  of  the  program  tracks  was  to  develop  Bio-SPICE  modeling 
framework,  which  stands  for  the  Simulation  Program  for  Intra-Cellular  Processes. 

During  the  period  of  January  2002  to  December  2003  CDRC  was  involved  in  the  DARPA 
BioCOMP  Program  as  a  subcontractor  to  the  Lawrence  Berkeley  Laboratory  (LBL)  team  led  by 
Dr.  Adam  Arkin.  Dr.  Philip  Colella  of  LBL  was  the  lead  scientist  in  charge  of  spatial  model 
development  in  the  LBL  Bio-SPICE  Team.  CFDRC’s  project  task  was  first  to  formulate  and 
explore  a  spatiotemporal  (multidimensional)  modeling  approaches  and  then  to  implement  and 
demonstrate  it  on  cell  and  organ  biology  problems. 

The  LBL  team  was  responsible  for  the  development  of  an  open  source  software  tool,  CHOMBO, 
for  spatiotemporal  modeling  of  cellular  biology  (http://seesar.lbl.gov/anag/chombo/index.html). 
CFDRC’s  task  was  to  formulate  fundamental  mathematical  models  and  their  numerical  solution 
procedures  for  solving  spatiotemporal  cell  biology  problems.  The  models  were  tested  using  the 
existing  CFDRC  software  framework,  CFD-ACE+  and  CFD  Micromesh.  Successful  models 
were  presented  to  the  LBL  team  for  further  implementation  in  the  Bio-SPICE  and  CHOMBO 
framework.  CFDRC  was  also  in  charge  of  the  development  of  3D  cell  membrane  model  that 
could  be  implemented  in  the  CHOMBO  code.  This  project,  completed  in  December  2003, 
provided  excellent  starting  point  for  the  follow  up  project  with  DARPA/IPTO  and  the 
AFRL/IFTC  Rome  NY. 
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Project  Objectives 

The  overall  objective  of  the  LBL-CFDRC  team  in  the  BioCOMP  Program  was  to  develop  a 
software  toolkit  and  modeling  framework  for  multidimensional,  spatiotemporal  comprehensive 
modeling  of  cell  biology.  The  objective  of  the  CFDRC  task  was  twofold: 

a)  Formulate  and  implement  multidimensional  cell  biology  models  in  CFD-ACE+  software 
and  to  perform  parametric  simulation  and  validation  tests  for  selected  spatiotemporal  cell 
biology  problems,  and 

b)  Develop  a  generalized  mathematical  model  of  a  cell  membrane  and  implement  it  in  a 
software  module  for  integration  in  the  spatiotemporal  cell  biology  model  developed  by  the 
LBL  team. 

Both  objectives  of  that  project  have  been  successfully  accomplished.  CFDRC  has  demonstrated 
exciting  simulation  results  for  several  challenging  cell  biology  problems  such  as  cellular 
chemosensing,  chemotaxis  and  cell  division.  C++  software  tools  for  3D  modeling  of  cell 
membrane  biology  has  been  developed,  documented,  tested,  and  delivered  in  open  source  form 
to  the  LBE  Bio-SPICE  team.  This  Appendix  presents  summary  of  results  of  that  project. 

Project  Approach 

During  this  project,  01/2002  -  12/2003,  CFDRC  has  formulated  fundamental  mathematical 
models  and  their  numerical  solution  procedures  for  solving  spatiotemporal  cell  biology 
problems.  The  models  were  tested  using  the  well-established  CFDRC  software  framework, 
CFD-ACE+.  Successful  models  were  presented  to  the  LBL  team  for  further  implementation  in 
the  Bio-SPICE  and  CHOMBO  framework.  Most  of  the  models  were  formulated  based  on 
models  published  literature.  Typically  such  models  were  formulated  as  OD  or  ID  problems. 
Our  approach  was  to  reformulate  these  published  models  as  2D  and  even  3D  transient  problems 
and  to  perform  parametric  simulations  to  assess  model  performance  in  the  multi  dimensional 
form.  C++  software  tools  for  3D  modeling  of  cell  membrane  biology  have  been  developed, 
documented,  tested,  and  delivered  in  open  source  form  to  the  LBL  Bio-SPICE  team.  The  results 
of  this  project  provided  solid  foundation  for  the  follow  up  project  directed  by  the  Air  Force 
Research  Laboratory  AFRL/IFTC. 

During  that  period  CFDRC  BioSPICE  team  has  been  involved  in: 
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Technical  discussion  with  BioSPICE  PI  and  team  members, 

Literature  review  on  biophysics,  biochemistry  and  modeling  of  Bacillus  Subtilis 
sporulation  and  chemotaxis  of  bacterial  and  immune  cells 
-  Automated  geometry  and  mesh  generation  for  realistic  cell  shapes. 

Formulation  of  mathematical  models  for  cell  membrane,  membrane  receptors,  cell 
transport  of  macromolecules  in  cytoplasm,  chemosensing,  chemotaxis  and  others 
Development  of  core  functions  for  0D/1D/2D/3D  cell  biophysics  processes 

The  following  sections  present  summary  of  results  obtained  during  the  course  of  the  CFDRC 
project  with  LBL  under  the  DARPA  BioCOMP  Program. 


Generation  of  Computational  Meshes  Spatiotemporal  Models  of  Cell  Biolosv 


Spatiotemporal  modeling  is  typically  performed  on  a  numerical  grid  (also  called  mesh)  forming 
set  3D  control  volumes  or  finite  elements.  The  generation  of  the  computational  grid  for  complex 
geometries,  such  as  3D  biological  cell,  is  a  challenging  problem.  In  this  project,  we  have 
explored  two  ways  to  generate  computational  grids  for  biological  cells:  image  based  mesh 
generation,  and  meshing  idealized  geometries  representative  of  biological  ells  e.g.  sphere  and 
rod.  Early  in  the  Bio-SPICE  project,  CFDRC  has  adapted  CFD-Micromesh  software  for  the 
application  in  automated  image  processing,  geometry  construction,  and  for  generation  of  2D  and 
quasi  3D  (extruded  2D)  meshes  for  simulation  of  cellular  biophysics  and  biochemistry.  Figure 
A.  1  presents  example  results  in  terms  of  original  cell  image,  unstructured  tetra  mesh,  adaptive 
binary  tree  mesh,  and  multi  block  structure  mesh  for  a  typical  shape  of  the  Bacillus  Subtilis  cell 
during  germination  process.  Note  the  cell  membrane  capturing  capability  and  formation  of 
septum. 


Figure  A.  1:  Bacillus  Subtilis  geometry  and  computational  mesh  during  germination  cell  cycle 
division.  Image  and  mesh  examples  obtained  with  CFD-Micromesh 


CFD-Micromesh  has  been  also  used  to  generate  mammalian  cells  with  embedded  nucleus. 
Geometric/Mesh  models  should  preserve  topological  information  about  membranes  and 
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organelles  so  that  separate  equations  (e.g.  for  signaling  or  sensing  receptors  or  membrane 
deformation)  can  be  solved.  Figure  A.2  shows  the  image-based  computational  mesh  with 
volume  conditions  for  the  nucleus,  cytoplasm,  and  plasma/nuclear  membranes. 


Figure  A.2:  CFD-Micromesh  generated  geometrical  image  and  computational  mesh  for  and 
immune  cell  with  a  complex  plasma  membrane  and  an  example  organelle  (nucleus). 

Model  of  Membrane  Geometry 

The  membrane  software  module  was  developed  by  CFDRC  and  was  delivered  in  mid  March 
2003  to  UCB/LBL  in  the  source  form.  The  membrane  module  was  written  in  C++  using  object- 
oriented  programming  with  a  very  powerful  database  for  geometry,  mesh,  solution  fields, 
properties,  and  boundary/volume  conditions.  It  combines  comprehensive  surface  constrained 
diffusion  models  with  simple  reaction  kinetics.  It  is  anticipated  that  the  kinetics  module  will  be 
replaced  with  a  much  more  advanced  reaction  kinetics  module  from  the  Bio-SPICE  LBL  team. 
This  team  will  take  the  lead  in  integrating  the  membrane  module  with  the  LBL  CHOMBO  main 
code.  CFDRC  has  performed  initial  tests  as  described  in  the  Results  section. 

It  has  been  decided  that  the  initial  version  of  the  membrane  module  will  be  driven  from  an  input 
text  file;  the  grid  will  be  read  from  a  file  with  the  CFDRC  MFC  format,  and  output  in  the  FAST 
format.  The  FAST  format  (developed  at  NASA  Ames)  is  an  “aerospace”  standard  format  for 
visualization  of  unstructured  grids  and  data  fields.  The  input  data  should  include:  mesh 
information,  boundary  and  initial  conditions,  and  solution  controls. 

The  membrane  module  code  delivered  to  LBL  consists  of  the  following  modules: 

In-core  database  with  all  variables,  properties,  mesh,  solution 
I/O  routines 

Main  driver  (not  needed  when  linked  with  LBL  solver) 

Interface  APIs  to  LBL  code 

Mesh  and  geometry  initialization  module 

Physical  data  initialization 

Mesh  movement  module 

Flux  calculassions  module  (diffusive,  advective) 

Simple  chemical  kinetics  source  terms  evaluation  module  (including  linearization) 

Matrix  assembly 
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Linear  equation  solvers 

Cellular  Chemosensing 

Cellular  ehemosensing  is  a  typical  spatiotemporal  cell  biology  problem.  It  involves  diffusion  of 
a  chemoattractant  in  3D  extracellular  matrix,  cell  membrane  binding  kinetics,  and  transport  of 
receptors  along  membrane  surface  toward  the  chemoattractant  gradient  -  all  coupled  to  an 
intracellular  signaling  network.  Figure  A.3  presents  an  example  chemosensing  problem 
proposed  by  Prof  Lauffenburger  team  at  MIT  (Narang  et  all  2001)  with  multi  step 
compartmentalized  chemical  kinetics  of  PIP3  signaling  in  cytoplasm  and  in  endoplasmic 
reticulum  (ER),  migration  of  membrane  receptors,  and  cell  polarization. 


Plasma 

Membrane 


Cytosol 


Endoplasmic 

Reticulum 


Figure  A.3:  Conceptual  model  of  ligand-receptor  binding  at  the  cell  membrane.  Schematics  of 
chemosensing  experiment,  and  biochemical  intracellular  IPS  signaling  model,  Narang  2001. 

In  this  simulation  the  chemoattractant  (ligand)  is  released  at  a  point  away  from  the  cell,  diffuses 
to  the  cell  membrane,  and  stimulates  membrane  bound  receptors.  The  receptors  respond  by 
diffusing  to  the  ligand  causing  receptor  crowding  in  that  zone.  We  solve  for  ligand  diffusion  in 
the  extracellular  space,  for  diffusion  and  crowding  of  membrane  bound  receptors,  and  for  the 
intracellular  signaling  pathway  as  sown  in  Figure  A.3  above.  The  cellular  membrane  has  been 
constructed  using  the  existing  membrane  model  with  membrane  equations  solved  only  for  the 
membrane  fixed  and  surface  diffusing  receptors. 

The  cell-membrane  model  includes:  3D  external  matrix  chemoattractant  diffusion,  binding  to 
membrane  receptors,  receptor  diffusion  along  membrane,  biochemical  reaction  on  membrane,  in 
cytosol,  and  in  external  matrix.  In  the  original  paper  by  Narang  and  Lauffenburger  the  model 
was  solved  as  one  dimensional  problem  (along  the  membrane  equator).  In  the  present  study  we 
model  as  2D  coupling  extracellular,  membrane,  and  intracellular  space. 
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Figure  A.4  and  Figure  A.5  present  example  simulation  results  for  the  cellular  chemosensing 
problem.  Some  of  the  variables  used  in  the  figures  below  are: 

•  rlO  :  Active  Receptors, 

•  71  :  Membrane  Phosphoinositides. 

•  TTs  :  Phosphoinositides  in  ER, 

•  i  :  Cytosolic  Inositol  and  Phosphates 

Figure  A.4  presents  transient  response  of  membrane  active  receptors,  rlQ,  and  membrane 
phosphoinositides,  ti,  to  an  increase  of  uniform  chemoattractant  concentration  for  two  levels  of 
stimulation:  stimulus  value  =  1  and  10.  Note  that  stronger  stimulation  solicits  more  intensive 
response  (red  lines).  Figure  A.5  presents  time  snapshots  of  membrane  phosphoinositides 
concentration  ti  during  the  stimulation. 


dark  line  stimulus  value  =1  ,  red  line  stimulus  value  =10 


Figure  A.4:  Transient  response  to  an  increase  of  uniform  chemoattractant  concentration  for  two 

values  of  chemoattractant  stimuli. 
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Figure  A.5:  Membrane  phosphoinositides  concentration  k,  as  a  function  of  cell  response  time. 
Top  left  initial  conditions,  bottom  right  final  equilibrium  state. 
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Cell  Cycle  Biochemistry  Model 

During  the  second  year  of  the  BioCOMP  Program  CFDRC  and  LBL  discussed  best  methods  for 
coupling  the  Bio-SPICE  biochemistry  problems  described  with  ODEs  and  multidimensional  cell 
biophysics  described  with  PDEs.  To  evaluate  computational  algorithms  needed  for  that  coupling 
CFDRC  has  stared  an  effort  to  link  simple  metabolic  pathways  (ODEs)  with  multidimensional 
cell  morphology  and  geometry  (PDEs).  We  have  chosen  Yeast  cell  cycle  model  as  a  test  case 
and  ACE+  Multiphysics  as  the  simulation  software. 

The  usual  form  for  the  cell  cycle  models  is  identification  of  relevant  CDKs,  catalytic  agents  and 
other  enzymes,  and  their  effects  on  various  stages  in  the  cell  cycle,  such  as  G1  phase,  S  phase 
and  so  on.  The  models  involve  rate  expressions  for  change  in  the  level  of  these  chemical  agents, 
in  the  form  of  ODEs,  which  are  interdependent  via  concentration  dependent  rate  constants.  We 
have  incorporated  a  robust,  coupled  ODE  solver  in  ACE+  for  the  modeling  of  cell  cycle, 
expressed  in  the  rate  expression  form.  As  an  example,  we  simulated  the  fission  yeast  cell  cycle 
developed  by  A.  Sveiczer,  John  Tyson  and  others  [Sveiczer  2001].  The  model  includes  eleven 
ordinary  differential  equations  and  several  rate  functions. 


We  simulated  the  ‘wild’  type  cell  cycle  as  well  as  several  mutant  cycles  involving  changes  in  the 
rate  constant  expressions  in  some  of  the  key  CDKs.  Figure  A.6  below  shows  the  calculated 
concentration-time  variation  of  several  of  the  CDKs  and  other  relevant  enzymes  during  the  cell 
cycle  of  the  wild  type  fission  yeast  cells.  These  compare  well  with  the  results  from  the  original 
reference  [Sveiczer  2001]. 

Table  below  presents  the  cell  cycle  model,  in  terms  of  the  rate  equations  (Sveiczer  2001): 


Table  A.  1:  Cell  Cycle  Model 

d 

kpf 
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Rate  functions 
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Figure  A.  6:  CFDRC  simulation  results  for  time-variation  of  various  chemical  species  in  the 

yeast  cell  cycle  model 

Multidimensional  Modeling  of  Bacterial  Sporulation  Biophysics 

The  common  accepted  model  to  understand  the  cell  division  process  is  focusing  on  the  behaviors 
of  Min  proteins  and  Z-rings  inside  cell  and  along  cell  membrane  [Meinhardt  2001,  Howard 
2001].  According  to  this  pole-to-pole  model,  the  proteins,  MinC,  MinD  and  MinE,  etc.  inside 
the  cell  oscillate  from  one  side  to  another.  The  MinE  protein  forms  in  the  center  of  the  cell  along 
the  membrane  at  the  beginning  of  oscillation.  Once  MinE  forms  it  will  chase  the  MinC/D 
proteins  staying  on  the  one  pole  of  cell.  Mean  time  ions  inside  cell  form  one  Z-ring,  Ftsz  ring, 
along  cell  membrane.  They  shrink  toward  the  center  point  of  cell.  MinC/D  will  be  pressed  to 
diffuse  into  other  pole  of  cell  slowly.  When  MinE  reaches  the  pole,  it  disappears  while  MinC/D 
located  on  the  other  side  of  cell  and  FtsZ  ring  divides  the  cell  into  two  identical  daughter  cells. 
This  process  will  continue  to  both  daughter  cells. 

Figure  A.  7a  below  shows  the  cell  structure,  mins  proteins  and  FtsZ  ring  and 

Figure  A.  7b  shows  the  formation  and  motion  of  FtsZ  ring  [Lutkenhaus  1997,  Erickson  2001]. 
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Figure  A.  7;  Image  of  bacterial  bell  early  in  the  sporulation  event  and  an  oscillatory  behavior  of 
Min  proteins  and  formation  FtsZ  ring,  [Lutkenhaus  1997,  Erickson  2001]. 


During  the  first  year  of  the  BioCOMP  Program  CFDRC  has  implemented  a  2D  model  of 
bacterial  sporulation.  We  followed  the  ID  approach  recently  proposed  for  FtsZ  ring  formation 
and  sporulation  for  E  Coli  [Howard  2001].  In  the  model  description  below  the  “Mins”  proteins 
bounded  by  membrane  are  denoted  by  D,  E  and  F  while  the  “Mins’  proteins  of  free  form  are 
denoted  by  d,  e,  and  f  The  letters  D  (d),  E  (e)  and  F  (/)  mean  MinD,  MinE  and  FtsZ 
respectively.  The  transportations  of  those  three  species  are  governed  by  the  following  reaction- 
diffusion  equations  in  multi  dimension  form: 
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Ot  OXi  OXi 

dd  2  N  j  ^ 

=  -  Pod^D  +cr£,)-//^J  +  — r^  — 

dt  dXi  dXi 


dt 

de 

~dt 


Pe^ 


D 


1  -l-  K£)^D^  1  -l-  k^E^ 


-  PeF  +^f  £  ^ 

dxi  dxf 


D  E  +  (Tp  d  „  de 

=  (Te-pEe - ^ - --p^eP  —  Ye  — 

\  +  Kp,eD^  \  +  KeE^  5xi  dXf 


(A.1) 

(A.2) 

(A.3) 

(A.4) 

(A.5) 

(A.6) 
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Where,  p  is  catalytic  production  rate,  a  is  constant  production  rate,  k  is  saturation  and  mutual 
interference,  |a  is  decay  rate,  F  is  diffusion  coefficient.  The  index  denotes  the  corresponding 
species.  The  distributions  of  proteins  will  be  given  by  solving  equations  A.l  to  A.6.  The  cell 
division  process  will  be  solved  by  geometry  deformation  approach.  Figure  A.  8  presents 
conceptual  representation  of  Bacillus  Subtilis  cell  division  process  and  a  2D  computational  mesh 
used  for  parametric  simulations. 


Figure  A.8:  E-Coli  cell  geometry  and  mesh. 

In  order  to  avoid  the  singularity  of  geometry  during  the  cell  deformation,  we  create  a  very  tiny 
flat  region  in  the  center  of  cell  where  FtsZ  ring  will  stay  and  shrink.  The  initial  protein  locations 
are  shown  below.  Figure  A.9  shows  the  oscillation  history  of  MinD  protein  and  cell 
deformation.  At  the  final  stage  the  septum  is  completely  closed. 


Figure  A.9:  Predicted  concentration  oscillation  history  of  MinD  protein  and  cell  deformation. 

Polarization  and  Oscillations  during  E  Coli  Cell  Division 

Bacterial  cell  division  is  driven  by  the  formation  of  a  contractile  FtsZ  protein  ring  at  the  cell 
division  site.  Proper  placement  of  the  FtsZ  ,  a  special  ring  formed  at  the  cell  midpoint  between 
segregating  daughter  chromosomes,  is  mediated  by  two  mechanisms,  1)  nucleoid  occlusion,  and 
2)  MinC,  Mind,  MinE  protein  system  spatiotemporal  oscillation.  In  E.  coli  during  each 
oscillation  period  MinD  accumulates  in  the  cell  membrane  in  a  “polar  zone”  at  one  end  of  the 
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cell.  This  polar  zone  then  shrinks  toward  the  end  of  the  cell  as  a  new  accumulation  forms  at  the 
opposite  pole.  MinC  follows  the  same  pattern  as  MinD,  and  the  two  proteins  form  complexes  on 
the  membrane.  In  contrast,  MinE  forms  a  ring  at  the  boundary  of  the  MinD  polar  zone  with 
some  MinE  dispersed  throughout  the  polar  zone.  The  MinE  ring  moves  toward  the  end  of  the 
cell  as  the  MinD  polar  zone  shrinks.  Figure  A.  10  presents  the  life  cycle  and  division  of  E.  coli. 


Figure  A.  10:  Schematic  and  a  microcopy  image  illustration  of  life  cycle  and  division  ofE.  coli. 

Three  models  of  the  Min  protein  oscillations  have  been  recently  proposed:  Howard  (2001), 
Meinhardt  (2001),  and  Huang  (2003).  All  of  these  models  successfully  generate  oscillations,  but 
none  can  be  meaningfully  compared  with  in  vivo  observations. 

Meinhardt  (2001)  model  used  six  ID  diffusion-reaction  equations  for  MinC,  D,  E  at  the 
membrane  and  Mine,  d,  e,  in  the  cytoplasm.  The  model  requires  newly  synthesized  protein  to 
form  both  the  MinD  polar  zones  and  MinE  ring,  because  the  proteins  disappear  (are  degraded) 
from  the  simulation  on  dissociation  from  the  membrane.  This  is  not  fully  confirmed  in 
experiments.  That  model  also  arbitrarily  assumes  that  MinE  attaches  preferentially  to  an 
intermediate  concentration  of  MinD. 

The  model  of  Howard  (2001)  produces  oscillations  only  if  MinE  is  driven  onto  the  membrane 
by  cytoplasmic  MinD,  despite  evidence  that  MinE  is  recruited  to  the  membrane-by-membrane- 
associated  MinD.  In  addition  the  model  predicts  oscillations  that  have  the  opposite  dependence 
of  frequency  on  MinD  concentration  than  is  observed,  and  MinD  forms  a  medial  band  that 
moves  toward  the  end  of  the  cell,  contrary  to  the  experimental  observation  that  MinD  forms 
directly  as  a  polar  zone. 

In  Hwang  (2003)  model  the  oscillations  are  driven  by  hydrolysis  of  ATP:  (1)  MinD  binds  ATP 
and  the  complex  binds  to  the  membrane;  (2)  MinE  binds  to  the  MinD: ATP,  induces  ATP 
hydrolysis,  and  all  constituents  are  released  from  the  membrane;  (3)  MinD:ADP  releases  ADP, 
completing  the  cycle.  This  is  the  only  model  using  cell  bioenergetics  of  ATP. 

The  equation  system  used  in  the  Meinhardt  de  Boer  (2001)  model  is  presented  in  table  below. 
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Note  that  the  membrane  bound  proteins  are:  MirrD,  E  and  FitsZ  and  the  free  intracellular  proteins 
are  Min^^,  e  and  f. 


Table  A. 2:  Equation  system  used  in  the  Meinhardt  de  Boer  (2001)  Model. 


FitsZ 

dE  F^  +  (7 

—  =  V(r^VF)  +  p,f- - ^  -  p,F  - 

Ot  1  +  KpF 

FitsZ 

^{=V(r,V/)  +  cr^  P,f 

MinD 

QD 

—  =  V(r^ VD)  +  +  cr^ )  -  PpP  -  jUp^pDE 

at 

MinD 

dd 

=  s;/(rjVd)  +  aj-pp,d(D^  +  (7p,)-ju^d 

at 

MinE 

8E  D  E^+(7p 

=  V{YpVE)  +  PEe  2  •  1  ^  biEE 

ut  1  +  Kj^^D  1  + 

MinE 

de  .  D  E^  +(7p 

=  V(r^Ve)  +  cr,  Ppe  •  p^e 

at  1  +  1  +  k^E 

Where  p  is  catalytic  production  rate,  a  is  constant  production  rate,  k  is  saturation  and  mutual 
interference,  and  p  is  decay  rate.  The  nonlinear  kinetics  of  binding  cytoplasmic  proteins  to  the 
membrane,  and  an  order  of  magnitude  larger  cytoplasmic  diffusion  compared  to  membrane 
diffusion  result  in  the  spatial  oscillations  of  Min  proteins. 

The  model  proposed  by  Howard  (2001)  uses  only  four  coupled  reaction  diffusion  equations 
describing  respectively  the  concentrations  of  proteins  MirrD  in  the  cytoplasm.  Mind  on  the 
membrane,  MinE  in  the  cytoplasm  and  Mine  on  the  membrane.  The  model  does  not  have 
explicit  conservation  equation  for  FitsZ.  Table  below  presents  the  equation  system  for  the 
Howard  model. 

Table  A.  3:  Equation  system  for  the  Howard  Model. 
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The  model  proposed  by  Huang  (2003)  is  also  using  only  MinD  and  MinE  proteins  but  it  includes 
the  ATP  hydrolysis  driven  membrane  binding  kinetic  mechanism.  The  specific  interactions  of 
the  Huang  model  are  shown  schematically  in  Figure  A.l  1.  The  oscillations  are  driven  by  a  cycle 
in  which  MinD: ATP  complex  first  associates  with  the  membrane,  preferentially  where  other 
MinD:ATP  is  located.  MinE  then  attaches  to  the  MinD:ATP,  activates  ATP  hydrolysis,  and 
MinE  and  MinD:ADP  reenter  the  cytoplasm  in  a  1:1  ratio.  The  equations  describing  the  time 
evolution  of  MinD  and  MinE  concentrations  in  a  cylindrical  cell  are  shown  in  Figure  A.l  1. 


4 


dD-.ADP 


dt 

dD-.ATP 


=  r^V"Z) :  ADP  -  k,D :  ADP  +  k^”de 


dt 


=T„VD:ATP-kfi\ADP-[k:;  +k,'”(d+dd)} 


—  =  r,V^E  +  k,,de  -  k"d  ■  E 

dt  E  2 

dd 

dde 
dt 


=  -k^dE  +  {k,  +k,{d  +  de)]  D:  ATP 


=  -k"de  +  k"d-E 


Figure  A.l  P-  Schematic  and  a  microcopy  image  illustration  of  life  cycle  and  division  ofE.  coli. 

Shown  are  equations  A.  7  thru  A.  10,  Huang  2003. 

During  this  project  we  have  used  both  CFD-ACE+  and  the  Membrane  Model  software  (later  also 
repeated  with  CoBi  tools)  to  test  all  three  models  and  we  have  reproduced  their  behavior.  The 
3D  model  was  used  to  assess  the  scalability  of  the  code  and  to  explore  potential  anomalous  and 
interesting  behavior  of  the  models  for  3D  bacterial  cell  geometries. 

Figure  A.  12  presents  example  simulations  results  for  E.  Coli  polar  oscillations  of  Min  proteins 
obtained  for  ID  and  for  3D  cell  geometries.  The  model  was  setup  as  a  transient  reaction- 
diffusion  problem  with  coupled  equation  systems  for  cytoplasmic  and  membrane  bound  Min 
proteins.  As  shown  in  Figure  A.  12  our  model  properly  predicted  unstable  oscillations  of  MinE 
and  MinE  proteins  and  focusing  of  FtsZ  protein  near  the  center  of  the  cell,  septum  formation  site. 
When  the  model  was  run  in  3D  configuration  with  a  rod  like  bacterium  we  have  noticed  an 
interesting  pattern  of  not  only  axial  instabilities  but  also  circumferential  instabilities  as  well. 
Figure  A.  12  b  presents  an  FtsZ  protein  concentration  on  the  bacterial  membrane. 

Notice  that  the  code  properly  predicts  maximum  concentration  in  the  middle  of  the  cells  but  it 
also  shows  that  FtsZ  is  initially  preferentially  focusing  at  a  point  on  the  membrane.  We  observe 
that  the  focus  point  is  circumferentially  moving  on  the  membrane  surface.  This  has  not  been 
reported  in  the  literature  before. 

The  results  obtained  with  the  Membrane  Model  software  have  later  been  successfully  reproduced 
with  the  CoBi  code.  CoBi  has  the  potential  to  predict  experimentally  observed  formation  of 
spiral  arrays  of  Min  proteins  on  the  inner  side  of  the  cell  membrane. 
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Figure  A.l  2:  ID  and  3D  simulation  results  for  the  bacterial  oscillations  of  Min  proteins  during 
bacterial  division.  Results  obtained  with  the  membrane  model  (later  reproduced  with  the  CoBi 

code). 

One  of  the  major  challenges  in  computational  spatial  cell  biology  is  the  solution-  and  time- 
dependent  geometry  of  the  dividing  cell.  We  decided  to  demonstrate  the  CoBi  code  on  modeling 
an  idealized  problem  involving  cell  division,  transient  diffusion  of  membrane  receptors  and 
simple  external  ligand  cell  receptor  binding  kinetics.  Figure  A.  13  presents  two  time  instants  of 
the  membrane  geometry  and  unstructured  deforming  mesh  just  prior  to  division  initiation  and  at 
the  cell  separation  instance.  Note  that  during  the  division  process  surface  receptors  diffuse 
between  two  parts  of  the  membrane.  Figure  A.  14  presents  selected  snapshots  of  the  membrane 
receptor  concentration  color  maps  during  the  cell  division. 


Figure  A.  13:  Bacterial  cell  division  simulation  with  the  membrane  model  of  the 
Bio-SPICE  cell  biology  framework. 


Figure  A.  14:  Transient  distributions  of  membrane  receptors  during  bacterial  cell  division. 
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Bacterial  Chemosensing  and  Chemotaxis 


During  this  project  CFDRC  and  LBL  selected  bacterial  chemotaxis  as  the  most  challenging 
example  of  spatiotemporal  cell  biology.  Figure  A.  15  schematically  illustrates  E.  Coli  chemotaxis 
signaling  pathways  and  their  spatial  intracellular  locations.  Mathematical  modeling  of  bacterial 
chemotaxis  would  need  to  integrate  several  biophysical  and  biochemical  phenomena  such  as: 
signaling  pathways,  energetics,  intracellular  diffusion,  molecular  motor  chemo-mechanics,  and 
flagella  propulsion  and  bacteria  dynamics  in  bio  films  with  spatial  chemoattractant  gradients. 
Following  is  the  description  of  the  salient  features  of  the  model  and  simulation  results  of  model 
components. 


Figure  A.  15:  Schematic  of  signaling  pathways  for  E  Coli  chemotaxis. 

[Eisenbach  2001,  Spiro  1997J. 

The  bacterial  motion  is  driven  by  the  rapid  rotation  of  helical  flagellar  filaments  that  protrude 
from  the  cell  wall.  Flagellum  rotation  is  driven  by  a  motor  located  in  the  cell  wall  at  the  aft  side 
of  the  cell.  The  motors  are  powered  by  the  chemo-mechanical  reactions  and  resultant  flux  of 
ions  (H+,  Na+)  across  the  plasma  membrane.  In  an  isotropic  environment  a  bacterium  swims 
about  in  a  random  walk  produced  by  alternating  episodes  of 

Runs,  counterclockwise  (CCW)  rotation  of  the  flagellum  and  motor,  and 
Tumbles,  clockwise  (CW)  flagellar  rotation. 

In  a  field  of  chemoattractant  gradient,  the  cell  uses  multiple  transmembrane  chemoreceptors  (e.g. 
methyl-accepting  chemotaxis  proteins  -MCPs)  to  bind  ligands  on  the  periplasmic  side.  Binding 
events  stimulate  cascades  of  cytoplasmic  signaling  reactions,  which  create  functionalized 
chemicals  that  extend  the  duration  of  runs  (duration  of  CCWs)  that  happen  to  carry  it  in 
favorable  directions. 

The  cytoplasmic  signaling  cascade  is  initiated  close  to  the  MCPs  where  several  proteins  (CheA, 
CheY,  CheZ)  are  activated  by  methylation  and  phosphorylation  events,  some  of  which  (CheY) 
move  to  the  aft  part  of  the  cell  to  communicate  with  the  flagellar  motors  via  a  phosphorelay 
sequence.  Figure  A.  15  (above)  illustrates  potential  mechanisms  of  chemotactic  signaling  in  E. 
coli.  Detailed  discussions  of  the  mechanism  and  mathematical  model  have  been  published 
[Spiro  1997]. 
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During  the  BioCOMP  Program  we  have  modeled  two  signaling  pathways  for  E  coli  proposed  by: 


{SPO}  Spiro,  Parkinson,  and  Othmer,  PNAS  July,  1997,  and 
{ASBT}  Almogy,  Stone,  and  Ben-Tal,  Biophysical  J.,  Dec.  2001. 

We  have  extended  the  ASBT  model  to  multidimensional  form  by  implementing  diffusive  flux  of 
CheY  and  CheYP  within  the  bacterial  cytoplasmic  space.  We  present  comparison  of  ODE  and 
PDE  models. 

SPO  Model 

The  signaling  mechanism  of  the  SPO  model  is  illustrated  in  Figure  A.  15  (above)  and 
schematically  in  Figure  A.  16.  The  corresponding  reaction  mechanism  is  summarized  in  the  table 
below: 


Table  A.4:  SPO  model  reaction  mechanism. 


Description 

Reaction  Labels 

Kinetics  and  rate  labels 

Methylation 

kla 

Michaelis-Menten  kinetics 

ri...,r4 

T2  +  R  T2R  — T3  +  R 

kib 

First-order  kinetics 

h-.-AA 

T2+R - =^T2+Bp 

Demethylation 

r-i...,r-4 

T3+Bp-!^T2+Bp 

Ligand  binding 

r5,r6,r7;r_5,r_6,r-7 

ks 

T2  _L  «  LT2 
k-5 

Autophosphorylation 

r8»---Al3 

^7^ 

00 

Phosphotransfer 

h 

T2p  +  B-!^T2Bp 
T2p+Y-!ii-^T2  +  Yp 

Dephosphorylation 

Bp-ki^B 

Yp  +  Z  >  Y  +  Z 

Here  T2,  T3,  T4  represent  methylated  states  of  the  Tar  chemoreceptors,  abbreviations  A,B,R,Z 
are  used  to  denote  Che  proteins,  and  subscript  p  denotes  phosphorylated  state.  In  addition  there 
are  Tar  receptor  ligation  reactions  on  the  periplasmic  side.  Figure  A.  16  presents  complete 
schematic  of  the  reaction  system.  L  is  the  external  stimulus  to  which  the  system  must  adapt. 
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Figure  A.  16:  Ligand-binding,  phosphorylation,  and  methylation  reactions  of  the  Tar-CheA- 
CheW  complex,  and  detail  of  the  phosphotransfer  reactions,  Spiro  1997. 

The  dependent  variables  to  be  solved  by  the  differential  equations  are: 

T2,T3,T4,  LT2,LT3,LT4,  T2P,T3P,T4P,  LT2P,LT3P,LT4P,  YP,  and  BP. 


Example  speeies  conservation  equations  are: 

1  rj-1 

-^=  -Lk5T2+k_5LT2  -  ^  +  kyT2p(Yc  -  Yp)+  ki,T2p(B,  -  Bp) 

Lig3.n(i  binding  /  tgIgs-sg  Phosphorylation  Phosphotransfer 

-kR - ^ - +Bk,F  (A.ll) 

KbiKa+T  ^ 

^ ^  Demethylation 
methylation 


We  write  the  equations  for  CheYp  and  CheBp  as 

dY. 


dt 

dt 


where  TP  is  the  level  of  total  phosphorylated  receptor  complex: 


Tp  =  fp  +  Lfp  +  fp  +  Lfp  +  fp  +  Lfp 


We  have  made  use  of  the  conserved  quantities 


=  ^2  +  Bf  +f+  Lf  +f+  Lf  +  fp  +  Lfp  +  fp  +  Lfp  +  fp  +  Lfp 


(A.12) 
(A.  13) 


(A.14) 

(A.15) 


Y^=Y  +  Yp  (5.17)  and  B^=B  +  Bp  (A.16) 

Note  that  CheY  and  CheYP  are  the  two  key  species  that  diffuse  through  the  cytoplasm  from  the 
sensing  (chemoreceptor)  side  of  the  bacteria  to  the  aft  (motored)  side. 
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ASBT  Model 


This  model  has  been  formulated  and  demonstrated  by  in: 

Gal  Almory,  Lewi  Stone  and  Nir  Ben-Tai,  “Multi-State  Regulation,  a  Key  to  Reliable  Adaptive 
Biochemical  Pathways”,  Biophysical  Journal  81,  pp3016-3028,  Dec.,  2001. 

The  Chemotaxis  signaling  mechanism  of  the  ASBT  model  is: 


Re  ceptor(active)  +  CheY  — — > 

CheYp  +KQceptor  {active) 

CheZ  +  CheYp  +  CheZ(active) — ^ 

2CheZ{active)  +  CheYp 

CheZ  (ae  five)  — 

CheZ 

CheZ  {active)  +  CheY^  — 

CheY  +  p  +  CheZ  {active) 

CheAs{free)  +  CheZ{free)  — ^ 

{CheAs  -  CheZ] 

{CheAs  -  CheZ]  ^ — 

CheAs{free)  +  CheZ{free) 

[CheAs  CheZ)  +  CheYp  — ^ 

CheY  +  p  +  {CheAs  -  CheZ] 

The  rates  equation  for  CheYp ,  {CheAs-CheZ}  and  CheZ(active)  are  as  follows: 


d{CheY^) 

dt 


A{CheY,^,^,  -  CheY^ )-b-CheY^[ {CheAs  -  CheZ}  +  k  ■  ] 


diCheAs  CheZ)  ^  ^^A  _  [CheZ,^^, 

at 


{CheAs  -  CheZ}^-  d  {CheAs  -  CheZ) 


dt 


e{CheZ^^,^^  +  s)CheY^{CheZf^^^  -  CheZ^^,J-  f  •  CheZ^ 


Where,  A,  b  etc.  are  model  constants. 

We  simulate  transport  of  CheYp  from  production  site  to  reduction  site  in  the  following  diffusion- 
reaction  process: 


0(C^  +  ^  (A.17) 

dt  dx  dx 

Where,  D  is  the  diffusivity,  R  is  reaction  rate  at  production  site  and  reduction  site  respectively. 
Figure  A.17  schematically  illustrates  spatially  polarized  locations  of  the  production  and 
consumption  of  CheYp  and  the  central  part  of  the  cell  where  CheYp  diffuses  through  cytoplasm. 
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Flagellar 


Figure  A.  17:  Bacterial  chemotaxis  signal  transduction  pathway  including  diffusive  transport. 


At  the  production  site,  the  governing  equation  for  CheYp  is 


5(CheYp) 

5t 


a  5(CheYp) 

=  ^  D - +  A(CheY,„,„  -  CheYp ) 


Ri=A(CheY,,t,i- CheYp) 
At  the  diffusion  link,  the  governing  equation  for  CheYp  is 

a(CheYp)  a^gfCheYp) 

di  dx  dy^ 

R2  =0 


At  the  reduction  site  the  governing  equation  for  CheYp  is 
5(CheYp)  a  d(CheYp) 

- ^  ^  -  b  •  CheYp  [  {CheAs  -  CheZ}  +  k  •  CheZ ] 

oi  OX  dx 


(A.  18) 
(A.19) 

(A.20) 

(A.21) 

(A.22) 


Let  vector  Y(3)  denote  the  current  (new)  concentration  of  CheYp  at  the  site  1,  2,  and  3.  Where 


the  number  1,  2,  and  3  means  production  site,  diffusion  link  and  reduction  site.  With  the  finite 
volume(or  finite  difference)  method,  we  have  the  matrix  for  CheYp : 


AY  =  Y“  +dtR 


where 


A  = 


1  +  dtD 
-dtD 


0 


-dtD 
l  +  2dtD 
-dtD 


0 

-dtD 
1  +  dtD 


Y°  is  the  is  the  concentration  of  previous  time  step  and  R  is  the  reaction  rate  defined  as: 


(A.23) 


Y“(l) 

'Ri' 

Y“  = 

Y“(2) 

R2 

Y”(3) 

R3 

(A.24) 


D  is  diffusivity  and  dt  is  the  time  interval. 

We  have  solved  the  above  matrix  equation  using  two  approaches: 

Full  1 -dimensional  model  using  mesh  do  represent  diffusive  fluxes,  and 

Three  Compartment  model  by  converting  diffusive  fluxes  form  PDE  form  to  ODE  form  for 

3  points  along  the  length  of  bacterium. 


87 


In  the  limiting  case  the  PDE  with  3  grid  cells  should  give  identical  results  as  the  multi 
compartment  ODE  model; 

Figure  A.  18  shows  the  time  history  of  CheYp  at  different  sites.  Figure  A.  19  shows  the 

comparison  between  ODE  solution  of  equation  (2)  and  diffusion-reaction  equations.  We  can  see 
that  the  production  site  and  reduction  site  clearly  make  difference.  The  diffusion-reaction 
equation  should  be  more  accurate  to  handle  the  transduction  pathways. 


Figure  A.  18:  Time  history  of  the  concentration  of  CheYp  at  three  sites  along  the  E  Coli  body. 


Figure  A.  19:  Comparison  between  ODE  solution  and  diffusion-reaction  solution. 
Figure  A.20  shows  the  time  snapshots  of  the  CheYp  diffusion  process. 
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Figure  A. 20:  Distribution  of  CheYp  within  thee  coli  cytoplasm  during  the  transduction  process. 

It  should  be  pointed  out  since  we  do  not  know  the  rate  constants  and  diffusivity  of  CheYp , 

therefore  the  solutions  shown  here  is  need  to  be  verified.  What  we  do  is  to  propose  a  concept  to 
simulate  the  process  in  transduction  pathway. 

Three  Dimensional  Integrated  Model  of  E  Coli  Chemotaxis 

CFDRC  under  the  DARPA  BioCOMP  Program  has  demonstrated  a  3D  bacterial  chemotaxis 
model  within  CFD-ACE+  framework.  The  model  combines  chemosensing,  intracellular 
signaling,  chemo-mechanics  of  molecular  motors,  switching  mechanism,  hydromechanics  of 
flagellar  beating,  and  the  dynamics  of  bacterial  body  due  to  propulsive  and  viscous  drag  forces 
and  moments. 

The  flagellar  filament,  a  thin,  helical  tube,  joins  the  flagellar  basal  body  via  a  hollow,  flexible 
hook,  which  acts  as  a  universal  joint  [Berg  2000].  The  hook  is  connected  to  the  rod,  a  straight 
tube  that  forms  the  driveshaft.  Three  rings  surround  the  rotating  rod:  MS  ring  in  the  cytoplasmic 
membrane,  P  ring  in  the  peptidoglycan  layer,  and  the  L  ring  in  the  outer  membrane. 

Each  of  the  six  or  so  motors  on  the  surface  of  an  E.  coli  is  a  cylindrical  structure  about  50  nm  in 
diameter,  is  powered  by  the  transmembrane  proton  motive  force,  and  rotates  at  more  than  100 
revolutions  per  second.  It  consists  of  over  40  different  kinds  of  proteins  of  which  MotA  and 
MotB  make  up  a  proton  channel  through  the  cytoplasmic  membrane.  The  motor-switch  complex 
embedded  in  the  MS  ring  consists  of  three  proteins  FliG,  FliM,  and  FliN  that  constitute  the 
directional  motor  switch  formed  in  a  characteristic  bell-shape  structure  known  as  the  cytoplasmic 
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C-ring.  The  switching  is  controlled  by  the  concentration  of  the  small  diffusible  cytosolic  protein 
CheY  P,  which  binds  to  FliM. 

Flagellar  motors  of  E.  coli  are  powered  by  protons  flowing  from  the  periplasm  into  the 
cytoplasm  down  an  electrochemical  gradient,  through  the  discrete  channels  (approx  8)  in  the 
clearance  between  stationary  Mot  and  rotating  Fli  complexes.  We  understand  the  structure  and 
function  of  the  flagellar  motor  but  not  very  much  about  the  details  of  how  it  works,  or  how  it 
accomplishes  abrupt  switching  from  forward  (CCW)  to  reverse  (CW)  rotation. 

In  molecular  and  cellular  biology  compact  elegant  mathematical  models  are  rare.  Theories  and 
models  are  described  in  a  narrative  language  based  on  laboratory  observations.  A  detailed 
multidimensional  mathematical  model  of  cellular  chemotaxis,  including  molecular  pathways  and 
biophysical  dynamics  does  not  exist.  Recent  advances  in  molecular  and  structural  cell  biology  as 
well  as  remarkable  interest  in  computational  biology,  resulted  in  series  of  exiting  reports 
demonstrating  successful  component  models  of  cell  biology.  Most  outstanding  progress  has 
been  achieved  in  computational  cellular  biochemistry  where  cell  metabolism  is  described  by 
large  systems  of  ODEs.  Some  of  them,  related  to  E.  Coli  chemotaxis  are  identified  below. 

To  our  knowledge  this  is  the  first  attempt  to  integrate  the  models  of  biochemical  signaling 
pathways,  metabolic  pathways,  molecular  motor  functionality,  flagella  rotation,  with  the  compete 
bacterial  body  hydrodynamics.  Moving  parts  of  the  molecular  proton  motor  are  submicroscopic 
and  immersed  in  a  viscous  medium  (water  or  lipid),  so  the  Reynolds  number  is  very  small  and 
viscous  dissipation  effects  should  be  included  along  with  electrostatic  interactions,  protein 
conformational  changes,  and  proton  binding  kinetic  events  as  it  cascades  down  in  the  channel 
from  the  periplasm  to  the  cytoplasm. 

Computational  modeling  of  bacterial  propulsion  has  been  investigated  half  a  century  ago  by  Sir. 
Geoffrey  Taylor  in  “The  Action  of  Waving  Cylindrical  Tails  in  Propelling  Microscopic 
Organisms”  Proc  Royal  Soc.  Vol  21 1,  pp  225-239,  Feb  1952.  Since  then  several  semianalytical 
models  have  been  presented,  not  much  better  than  original  Taylor’s  work.  Recently  Japanese 
bio-nano  consortium,  ERATO,  led  by  Prof  Keiichi  Nanda  have  presented  first  2D  axisymmetric 
model  of  motile  bacterium  with  a  single  flagellum  using  full  Navier-Stokes  equestrians  and 
moving  deforming  grid  (personal  communication  between  Dr.  A.  Przekwas  CFDRC  and 
Professor  Keichi  Namba  2005,  also  http://www.nanonet.go.ip/english/mailmag/2004/011a.html. 

During  the  second  year  of  the  BioCOMP  Program  CFDRC  has  integrated  a  complete  E  Coli 
chemotaxis  model  coupling  signaling  pathways,  flagella  motors,  beating  and  rotating  flagella, 
and  bacterial  swimming  in  a  3D  fluid.  The  present  model  involves  overset  mesh  for  bacterial 
body,  molecular  motor  and  switching  sub-models,  flagella  rotation  and  large  displacement 
during  CCW-CW  switch. 

We  have  performed  several  simulations,  created  animations,  and  observed  bacterial  runs  and 
tumbles.  Figure  A.21  presents  the  model  definition  and  example  simulation  results  for  a 
bacterium  with  a  single  flagellum  and  with  four  flagella.  Figure  A.22  presents  computational 
simulation  results  for  bacterial  micro-hydrodynamics,  pressure  maps  and  vorticity  isosurfaces  at 
selected  time  instants  during  the  bacterial  run-tumble  motion. 
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Figure  A.21:  Computational  model  of  Motile  Bacterium  with  rotating  flagella  using  full  Navier 
Stokes  equations  micro-hydrodynamic  model  and  an  overset  mesh. 


Figure  A.22:  Computational  simulation  results  for  bacterial  micro-hydrodynamics.  Pressure 
maps  and  vorticity  isosurfaces  at  selected  time  instants  during  run-tumble  motion. 
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Multivalent  Binding 


Interaction  of  cells  with  other  cells,  bacterial  membrane  attachment,  virus  endocytosis,  and  cell 
interaction  with  functionalized  substratum  involves  multivalent  binding  (MVB).  Figure  A.23 
schematically  illustrates  the  concept  of  multivalent  ligands  (6  and  12  valent)  binding  to  single 
and  bivalnet  receptors  on  cell  or  substrate  wall. 


Figure  A.23:  A  6  valent  ligand  binding  through  3  sites  to  single  valence  cell  membrane  receptor, 
and  single  multivalent  ligand  bound  to  bivalent  cell  membrane  receptor  (e.g.  Ab). 


Multivalent  Ligand-Receptor  (L-R)  binding  can  be  mathematically  represented  as  n-th  order 
chemical  kinetics  (Laufenburger  1993;  Perelson,  1994).  Figure  A.24  presents  the  MVB  modeling 
concept  for  a  simple  bivalent  ligand  and  for  a  cell.  Note  that  the  binding  becomes  stronger  with 
more  sites  attached.  Biological;  cells  (bacteria,  immune  cells)  are  using  this  mechanism  for 
strong  attachment  to  substrates  in  the  flowing  media. 


Figure  A.24:  Chemical  kinetics  model  of  a  binding  of  a  globular  multivalent  ligand  (or  a  cell)  to 

trivalent  surface  receptors 

To  date  MVB  events  have  been  modeled  in  ODE  framework  using  n-order  kinetics.  Until  this 
project  only  monovalent  ligand-receptor  (L-R)  binding  model  was  available  in  CFD-ACE+.  In 
this  project  we  have  impenitent  MVT  in  both  ODE  and  PDE  framework.  It  will  be  demonstrated 
on  cell-cell  and  cell-surface  attachment. 

When  a  multivalent  ligand  of  valence  ‘n’  binds  to  a  cell  surface,  it  can  attach  to  the  surface  at  i  = 
l,...n  sites.  This  different  configuration  is  called  as  ‘molecule  bound  state’.  Let  Xn(i)  be  the 
concentration  of  the  bound  state  in  which  a  ligand  of  valence  n  is  bound  to  a  cell  surface  receptor 
at  the  i*'’  binding  site.  The  general  equation  that  governs  the  binding  reaction  rate  is  given  as 
(Sulzer,  1996;  Perelson,  1998): 

V„(i)  =  (n-/+l)^,RV„(/-l)+(i+l)A:^,.,,,X„(/+^  \<i<n-\ 

Kin)  =  k„RXnin-\)-nk_„Xn{n),  i  =  n 
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Where  kj ,  k_j  are  the  forward  and  backward  (or  association  and  dissociation)  rates,  and  R  is  the 
number  of  receptor  sites  on  the  cell  surface.  Both  the  total  number  of  receptor  sites,  Rt,  and  total 
number  of  ligands,  Xnj,  are  conserved.  The  constraint  restrictions  on  the  above  equation  are: 

R,  =R  +  f^iX„ii) 

i=l 

^nT  =  X„(0)  +  «fiX  ^niO 

i  =  l 

where  B  is  number  of  cells  per  unit  volume. 

If  we  define  the  various  concentrations  as  : 

Volume  concentration  :  mol/L, 

Cell  concentration  :  (mL)  , 

Surface  concentration  :  cm’ , 
then 

a  =  10^A/NA  cm^  mol. 


2 

where  A  is  the  cell  mean  surface  area  in  cm ,  Na  the  Avogadro’s  number, 
=  6.023  X 10  /mol .  Using  these  definitions,  Rj  is  receptor  density  per  cm  .  For  a  cell  with 
radius  of  3.5  pm,  the  surface  area  A  is  1.5*10’^  cm^,  and  the  value  of  a  «  2.5  x  10”^^ . 


We  further  define  cross-linking  affinity,  Kj  = 
dimensionless  cross-linking  affinity  =  K^R 


i  >  1 ,  binding  affinity  K  =  Kj  = 


and 


We  simulated  a  test  case  given  by  Sulzer  (1996)  involving  binding  between  an  8 -valent  ligand 
and  receptors,  both  for  steady  state  and  transient  cases.  We  performed  simulations  for  different 
values  of  the  key  parameters  k^and  B,  where  kx  refer  to  the  cross-linking  affinity  of  the  i*’’ 
valence.  For  the  present  calculations,  all  kx  for  all  valences  were  assumed  to  be  the  same.  Figure 
A.25  shows  the  calculated  binding  states,  where  i  is  valence  (on  the  abscissa),  and  Xg(i)  is  the 
binding  surface  density  for  valence  I  (on  the  ordinate).  Results  for  a  range  of  values  of  kx  are 
plotted  in  Figure  A.25;  the  remaining  parameters  in  these  calculations  are  outlined  in  the  figure 
caption.  We  can  see  that  the  binding  density  strongly  depends  on  the  parameter  k^.  This  is 
natural,  since  k^  is  the  ratio  of  forward  and  backward  association  and  dissociation  rates. 


Figure  A.25:  Distribution  of  bound  states  in  the  excess  receptor  regime  with  different  . 
Parameters:  K  =  10*M“‘,St  ^nXn^  =10‘°cm“^a  =  10“^'’cm2mol,B  =  10’/mL.. 
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Figure  A.26a  and  Figure  A.26b  show  the  transient  solutions  for  two  of  the  cases  shown  in  Figure 
A.25.  The  transient  results  are  thus  consistent  with  the  steady  results  shown  above  for  large 
times.  For  small  cross-linking  affinity  coefficients,  (e.g.  0.1),  concentrations  of  the  single  bond 
are  highest,  while  the  8-bond  concentrations  are  highest  when  the  affinity  coefficient  ratio  is 
large.  For  the  middle-range  values  of  the  affinity  coefficients,  the  peak  concentrations  are  at  the 
mid-range  of  the  bond  numbers,  in  the  present  case  at  4-valence,  with  the  single-  and  8-bond 
concentrations  both  very  low. 


Figure  A.  26:  Time  history  of  binding  density  for  each  state  at  two  different  kx  values.  Circles  are 

steady  state  solutions. 

A  methodology  for  modeling  of  multivalent  binding  (MVB)  described  above  was  implemented 
in  a  biochemistry  module  (in  CFD-ACE+).  This  methodology  is  based  on  solutions  of  ordinary 
differential  equations  (ODEs)  that  describe  the  rates  of  change  of  concentrations  of  different 
chemical  species.  We  have  also  coupled  the  MVB  calculations  with  the  particle  transport  model 
to  calculate  binding  forces  on  cells  as  they  attach  to  surfaces  and  cells  using  MVB  of  surface 
receptors. 
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APPENDIX  B 


Input  Data  for  Selected  Computational  Cell  and  Tissue  Biology 
Test  Cases  Performed  with  the  CoBi  Code  and 
Presented  in  Chapter  5  of  this  Report. 


This  Appendix  presents  annotated  input  fdes  ,  *.SIM,  for  test  cases  computed  with  the  CoBi 
code  and  described  and  discussed  in  Section  5  of  this  report.  Whenever  possible,  schematics  of 
the  test  cases  are  also  provided  for  ready  reference  and  easier  interpretation  of  the  input  fdes. 
After  completed  simulations,  CoBi  creates  an  output  file  for  visualization  of  the  solution  results. 
The  output  fde  is  created  in  the  vtk  format  (www.vtk.org/pdf/fde-formats.pdf).  The  vtk  fde  can 
be  viewed  using  Visit  visualization  tools  developed  by  the  Lawrence  Livermore  National 
Laboratory  (www.llnl.gov/visit/),  or  other  visualization  software. 

The  line  numbers,  01,  02,  ...  shown  in  the  input  fdes  for  examples  below  are  not  needed  in  the 
*.SIM  fde  used  by  the  CoBi  code.  Here  they  are  used  only  to  annotate  and  explain  individual 
commands  of  the  *.SIM  fde. 


Diffusion  Problem  on  Planar  Square  Thin  Film  - 

01  thickness  0.01 

02  iteration  10 

03  module  share 

04  VC  Membrane  const_dens  1 

05  module  species  fl 

06  sweeps  500 

07  relaxation  0. 1 

08  VC  Membrane  const_diff  0. 1 

09  be  Left  fix_value  0.0 

10  be  Right  fix  value  1.0 

11  be  Top  fix  flux  0.0 

12  be  Bottom  fix  flux  0.0 

13  be  Membrane_outside  fix_flux  0 

14  be  Membrane_inside  fix_flux  0 

In  the  input  fde  keyword  thickness  is  the  membrane  thickness.  Keyword  iteration  is  the  number 
of  iterations  for  the  steady  state  case  or  each  time  step.  Keyword  transient  is  used  to  set  the 
number  of  time  steps,  size  of  a  time  step  and  output  frequency.  Keyword  module  is  used  to 
specify  input  data  for  the  module.  Keyword  share  means  that  input  for  shared  data,  such  as 
density. 
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Keyword  vc  is  used  to  specify  a  volume  condition  followed  by  a  volume  name,  constant  density, 
and  the  value.  Keyword  species  is  used  for  input  of  species  module.  Keyword  sweeps  is  used  to 
specify  the  maximum  number  of  sweeps  for  the  species.  Keyword  relaxation  is  used  to  input 
linear  relaxation  between  a  previous  value  and  a  new  value.  The  next  line  vc  sets  the  constant 
diffiisivity  for  volume  Membrane.  Keyword  be  means  the  input  for  boundary  condition  followed 
by  boundary  name,  boundary  condition  type  (fixed  value/fixed  flux). 

Line  01  sets  membrane  thickness  to  0.01; 

Line  02  sets  number  of  sub-iteration  to  10; 

Line  03  declares  the  start  of  share  module; 

Line  04  sets  density  to  1  for  membrane; 

Line  05  declares  the  start  of  species  fl  equation; 

Line  06  sets  number  of  iterations  for  linear  equation  solver; 

Line  07  sets  relaxation  factor  for  fl  to  0.1; 

Line  08  sets  the  volume  condition  (vc)  of  species  diffiisivity  to  0.1  for  membrane; 

Line  09-10  set  fix  concentration  0  on  the  Left  and  Ion  the  Right  Boundaries; 

Line  1 1-12  set  zero  gradient  concentration  on  bottom/bottom  surfaces; 

Line  13-14  set  no  flux  boundary  condition  to  membrane  inside/outside  surfaces; 


Diffusion  on  a  Spherical  Cell  Membrane 

01  thickness  0.01 

02  iteration  10 

03  module  share 

04  vc  Membrane  const_dens  1 

05  vc  Top  const_dens  1 

06  vc  Bottom  const_dens  1 

07  module  species  cl 
08  vc  Membrane  const_diff  10 

09  vc  Top  const_diff  10 

10  vc  Bottom  const  diff  10 

11  be  Membrane_outside  fix_flux  0.0 

12  be  Top_outside  fix_value  1 .0 

13  be  Top_inside  fix_value  1.0 

14  be  Bottom  outside  fix_value  0.0 

15  be  Bottom  inside  fix  value  0.0 

16  module  species  c2 

17  vc  Membrane  const_diff  1 

18  vc  Top  const_diff  1 

19  vc  Bottom  const  diff  1 

20  be  Membrane_outside  fix_flux  0.0 

21  be  Top  outside  fix  value  0.0 

22  be  Top_inside  fix_value  0.0 

23  be  Bottom  outside  fix  value  1.0 


96 


24  be  Bottom_inside  fix_value  1.0 

In  the  file,  keyword  “thickness  ”  is  the  membrane  thickness.  Keyword  '‘iteration  ”  is  the  number 
of  estimated  iterations  to  solve  a  steady  state  case  or  for  each  time  step  for  the  transient  problem. 
Keyword  "transient”  is  used  to  set  the  number  of  time  steps,  size  of  a  time  step  and  output 
frequency.  Keyword  "module  ”  is  used  to  specify  input  data  for  the  module.  Keyword  "share  ” 
means  input  for  shared  data,  such  as  density. 

Keyword  “vc”  is  used  to  specify  a  volume  condition  followed  by  volume  name,  constant  density 
and  the  value. 

Keyword  "species"  is  used  for  input  of  species  module. 

Keyword  "sweeps"  is  used  to  specify  the  maximum  number  of  sweeps  for  the  species.  Keyword 
"relaxation"  is  used  to  input  linear  relaxation  between  the  previous  value  and  a  new  value. 

The  next  line  "vc  ”  sets  the  constant  diffusivity  for  volume  Membrane. 

Keyword  "be"  means  the  input  for  boundary  condition  followed  by  boundary  name,  boundary 
condition  type  (fixed  value/fixed  flux). 

Line  01  sets  membrane  thickness  to  0.01; 

Line  02  sets  number  of  sub-iteration  to  10; 

Line  03  declares  the  start  of  share  module; 

Line  04-06  set  density  to  1  for  membrane,  top  and  bottom  of  membrane  regions; 

Line  07  declares  the  start  of  species  cl  equation; 

Line  08-10  set  diffusivity  to  1  for  the  membrane,  top  and  bottom  of  membrane  regions; 

Line  1 1  sets  no  flux  boundary  condition  to  membrane  outside  surface; 

Line  12-13  set  fix  concentration  1  to  top  inside  and  top  outside  surfaces; 

Line  14-15  set  fix  concentration  0  to  bottom  inside  and  bottom  outside  surface; 

Line  16-24  repeat  procedure  of  Line  07-15  for  species  c2,  but  the  concentration  values  are  set  to 
the  values  opposite  to  cl. 


Diffusion-Reaction  of  a  Flat  Membrane 


01 

thickness  1.00 

02 

iteration  10 

03 

transient  100  0.01  10 

04 

module  share 

05 

vc  Membrane  const_dens  1 

06 

module  species  test_fl 

07 

sweeps  500 

08 

relaxation  0. 1 

09 

vc  Membrane  const_diff  0.1 

10 

be  Left  fix_value  0.0 

11 

be  Right  fix  value  1.0 

97 


12  be  Top  fix  flux  0.0 

13  be  Bottom  fix  flux  0.0 

14  be  Membrane_outside  fix_flux  0 

15  be  Membrane_inside  fix  flux  0 

16  module  speeies  test_f2 

17  ve  Membrane  eonst_diff  0 . 1 

18  be  Left  fix_value  1.0 

19  be  Right  fix  value  0.0 

20  be  Top  fix  flux  0.0 

21  be  Bottom  fix_flux  0.0 

22  be  Membrane_outside  fix_flux  0 

23  be  Membrane_inside  fix_flux  0 

Line  01  sets  membrane  thiekness  to  0.01; 

Line  02  sets  number  of  sub-iteration  to  10; 

Line  03  sets  number  of  time  steps,  time  step  size  and  output  frequeney; 

Line  04  deelares  the  start  of  share  module; 

Line  05  sets  density  to  1  for  membrane; 

Line  06  declares  the  start  of  species  _fl  equation; 

Line  07  sets  number  of  sweeps  for  the  linear  solver; 

Line  08  sets  relaxation  for  _fl  to  0.1; 

Line  09  sets  species  diffusivity  to  0. 1  for  the  membrane; 

Line  10-11  set  fix  concentration  to  0  for  left  boundary  and  to  1  for  right  boundary; 

Line  12-13  set  no  flux  boundary  condition  to  top  and  bottom  boundaries; 

Line  14-15  set  no  flux  boundary  condition  to  outside  and  inside  membrane  boundaries; 

Line  16-23  repeat  procedure  of  Line06-15  for  species  _f2,  but  the  concentration  values  are  set  to 
the  values  opposite  to  _fl . 


Multicellular  Perfused  Tissue  EGF  Signaling 


001 

002 

003 

004 

005 

006 

007 

008 

009 

010 


iteration  3 
transient  100  800  1 

DTP 

membrane  MEM  2e-9  Cell_Gap 


module  share 
VC  Cell  const_dens  1000 
VC  Gap  const_dens  1000 
VC  MEM  Cell  Gap  const  dens  1000 
VC  MEM_Cell_Gap  reaction  SignalModell  Mem  euler 
VC  Cell  reaction  SignalModell_Vol  euler 
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oil  module  flow  f 
012  sweeps  30  30 

013  diagrelaxation  .5  .5  .5 

014  relaxation  0  0  0  0.2 

015  VC  Cell  const_diff  0. 

016  VC  Gap  const  diff  1.57e-3 

017  VC  MEM_Cell_Gap  const  diff  0. 

018  be  Gap_Inlet  fix_velocity  0.001  0  0 

019  be  Gap  Outlet  fix_pressure  0.0 

020  be  Cell  symmetry 

021  be  Gap  symmetry 

022  be  MEM_Cell_Gap_Gap  fix  veloeity  0  0  0 
023  be  MEM_side_boundary  symmetry 

024  module  species  L 

025  sweeps  80 

026  VC  Cell  const_diff  0. 

027  VC  Gap  const_diff  l.e-16 

028  VC  MEM_Cell_Gap  const  diff  l.e-16 

029  be  Gap  inlet  fix_value  l.e-9 

030  be  Gap  Outlet  fix  value  0. 

03 1  be  Cell  fix  flux  0 

032  be  Gap  fix_flux  0 

033  be  MEM_Cell_Gap_Cell  cg_wall 

034  be  MEM_Cell_Gap_Gap  cg  wall 

035  be  MEM_side_boundary  fix_flux  0 

036  module  species  Li 
037  sweeps  80 

038  VC  Cell  const_diff  l.e-16 

039  VC  Gap  const_diff  0. 

040  VC  MEM_Cell_Gap  const  diff  l.e-16 

041  be  Gap  inlet  fix_value  0. 

042  be  Gap_Outlet  fix_value  0. 

043  be  Cell  fix_flux  0 

044  be  Gap  fix_flux  0 

045  be  MEM_Cell_Gap_Cell  cg  wall 

046  be  MEM_Cell_Gap_Gap  cg  wall 

047  be  MEM_side_boundary  flx_flux  0 

048  module  species  Rs 

049  sweeps  80 

050  VC  Cell  const_diff  0. 

051  VC  Gap  const  diffO. 

052  VC  MEM_Cell_Gap  const  diff  l.e-16 

053  be  Gap  inlet  fix_value  0. 
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054  be  Gap  Outlet  fix  value  0. 

055  be  Cell  fix  flux  0 

056  be  Gap  fix_flux  0 

057  be  MEM_Cell_Gap_Cell  eg_wall 

058  be  MEM_Cell_Gap_Gap  eg  wall 

059  be  MEM_side_boundary  fix_flux  0 

060  init  MEM_Cell_Gap  fix_value_bbox  0.05  -2-2  -2222 

061  module  speeies  Cs 

062  sweeps  80 

063  ve  Cell  eonst_diff  0. 

064  ve  Gap  eonst  diff  0. 

065  ve  MEM_Cell_Gap  eonst  diff  1  .e- 1 6 
066  be  Gap  inlet  fix_value  0. 

067  be  Gap_Outlet  fix_value  0. 

068  be  Cell  fix_flux  0 

069  be  Gap  fix_flux  0 

070  be  MEM_Cell_Gap_Cell  eg_wall 

071  be  MEM_Cell_Gap_Gap  eg  wall 

072  be  MEM_side_boundary  fix_flux  0 

073  module  speeies  Ri 

074  sweeps  80 

075  ve  Cell  eonst_diff  1  .e- 1 6 

076  ve  Gap  eonst_diff  0. 

077  ve  MEM_Cell_Gap  eonst  diff  1  .e- 1 6 
078  be  Gap_Inlet  fix_value  0. 

079  be  Gap_Outlet  fix_value  0. 

080  be  Cell  fix_flux  0 

081  be  Gap  fix  flux  0 

082  be  MEM_Cell_Gap_Cell  eg_wall 

083  be  MEM_Cell_Gap_Gap  eg  wall 

084  be  MEM_side_boundary  fix_flux  0 

085  init  MEM_Cell_Gap  fix_value_bbox  1 .68e-7  -2-2-2222 

086  module  speeies  Ci 
087  sweeps  80 

088  ve  Cell  eonst_diff  1  .e-16 

089  ve  Gap  eonst_diff  0. 

090  ve  MEM_Cell_Gap  eonst  diff  1  .e- 1 6 
091  be  Gap_Inlet  fix_value  0. 

092  be  Gap  Outlet  fix_value  0. 

093  be  Cell  fix  flux  0 

094  be  Gap  fix  flux  0 

095  be  MEMCellGapCell  egwall 

096  be  MEM_Cell_Gap_Gap  eg  wall 

097  be  MEM_side_boundary  fix_flux  0 
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Line  01  sets  number  of  sub-iteration  to  3; 

Line  02  set  number  of  time  steps  to  100,  time  step  to  800  and  output  frequency  to  1  time  step; 
Line  03  sets  input  file  as  DTP  format  (DTP  is  CPDRC  database  format); 

Line  04  creates  membranes  between  cell  and  gaps  and  sets  membrane  thickness  to  2e-9  meter; 
Line  05  declares  the  start  of  share  module; 

Line  06-08  set  density  to  1000  for  cells,  gaps  between  cells  and  membrane 

Line  09  sets  reaction  to  membrane  signaling  model  for  membrane  and  implicit  Euler  method  for 

ODE  solver; 

Eine  10  sets  reaction  to  cell  signaling  model  for  cells  and  implicit  Euler  method  for  ODE  solver; 
Line  1 1  declares  the  start  of  flow  module; 

Line  12  sets  number  of  sweeps  to  30  for  both  momentum  and  continuity  equations; 

Line  13  sets  inertia  relaxations  to  0.5  for  momentum  equations  in  x,  y,  z  directions; 

Line  14  sets  relaxations  to  0.2  for  pressure  equation  and  0  to  momentum  equations; 

Line  15-17  set  viscosity  to  1.57e-3  for  gaps  in  between  cells  and  0  for  inside  cells  and 
membranes(means  no  flow); 

Line  18  sets  fixed  velocity  to  0.001  in  x  direction  for  boundary  Gap_Inlet; 

Line  19  sets  fixed  pressure  boundary  condition  to  0  for  boundary  Gap_Outlet; 

Line  20-21  set  top  and  bottom  of  cells  to  symmetrical  boundaries; 

Line  22  set  non-slip  wall  condition  for  outside  wall  of  membranes; 

Line  23  sets  symmetry  boundary  for  rest  of  boundaries; 

Line  24  declares  the  start  of  species  Ligand  equation; 

Line  25  sets  number  of  sweeps  to  80  for  linear  equation  solver; 

Line  26-28  set  species  diffusivity  to  0  for  inside  cell  and  le-16  for  gaps  in  between  cells  and 
membrane; 

Line  29-30  set  Ligand  concentration  to  le-9  at  inlet  and  0  at  outlet; 

Line  31-32  set  no  flux  boundaries  at  top  and  bottom  of  cells; 

Line  33-34  set  permeable  wall  for  membrane  inside  and  outside  wall; 

Line  35  sets  no  flux  boundary  for  rest  of  boundaries; 

Line  36-97  repeat  procedure  of  Line24-35  for  species  Li,  Rs,  Cs,  Ri,  Ci,  but  the  diffusivities  and 
boundary  conditions  are  set  to  the  different  values. 

Morphogenesis  of  Yeast  Cell 

iteration  5 

transient  11500  0.01  50 
DTP 

movingGrid 


reaction  Y  east_Morphogenesis_Model 
method  euler 
abs_tolerance  l.e-1 
rel  tolerance  l.e-1 
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t  =  mod(time,150) 
ks  =  t<115?0.2:0.1 
kd  =  0.05 

k31  =t<115  ?6.  :  0. 
k32  =  t<115?8.  :0. 
k4  =5. 

ksu  =  t<  115  7  0.05:  0.025 
kdu  =  0.05 

vl  =  (k31+k32*Fu*Ffactin*Ffactin)*Fgactin 

v2  =  k4*Ffactin 

Ffactin'  =  -kd*Ffactin  +  vl  -  v2 

Fgactin'  =  ks  -  kd*Fgactin  -  vl  +  v2 

Fu'  =  ksu  -  kdu*Fu 

reaction  Y  east_Morphogenesis_Model_Mid 
method  euler 
abs_tolerance  5.e-l 
rel_tolerance  5.e-l 

t  =  mod(time,150) 
ks  =  t<  115  7  0.2:  0.1 
kd  =  0.05 

k31  =t<115  7  6.  :t>125  7  50.  :0. 
k32  =  t<115  7  8.  :0. 
k4  =5. 

ksu  =  t<  115  7  0.05  :  0.025 
kdu  =  0.05 

vl  =  (k31+k32*Fu*Ffactin*Ffactin)*Fgactin 
v2  =  k4*Ffactin 

Ffactin'  =  -kd*Ffactin  +  vl  -  v2 
Fgactin'  =  ks  -  kd*Fgactin  -  vl  +  v2 
Fu'  =  ksu  -  kdu*Fu 

module  share 
VC  VC  const_dens  1 

VC  VC  reaction  Yeast_Morphogenesis_Model 

module  species  Ffactin 
sweeps  30 
VC  VC  const_diff  1. 
be  BC  fix_flux  0 
init  VC  fix  value  0.0 
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module  species  Fgactin 
sweeps  30 

VC  VC  const_diff  120 
be  BC  fix_flux  0 
init  VC  fix  value  0.0 

module  species  Fu 
sweeps  30 

VC  VC  const_diff  0.5 
be  BC  fix  flux  0 
init  VC  fix  value  0.0 

Rate  Equations 


1  4-1  yps  r 

YGL^GP  ^MAX,GL^GP  y 


^GY^GP  ~ 


1 


-Cr. 


GY->-GP  '‘'MAX,GY^GP  „  /  \2  GY 

PS  [cs  y 


3 .  ^GP^PY  ~  A 


-Cr 


GP^PY  '^MAX,GP^PY  r  \2  GP 

, ,  w.  ,  f 

PS  [  l/RS  y 


4. 


^GP^GY  ~  ^MAXfiP^GY  ' 


1 


-c. 


IPS 


/  ,  \2  ^GP 

^1/CV^ 


V  yes  j 


5.  ^py^aC 


-c. 


PY^AC  '^MAX,PY^AC  \j  J^p  ^  RS 

l/AF  ^  l/RS 


6.  ‘PpY^LA  ~ 


RS 


-C, 


PY^LA  '''MAX,PY^LA  r>c 

Ko„  +  Kij 


■  ^AC^CO^  ~  \ 


-C. 


AC^C02  ''‘MAX  ,AC^C0^  l/RS 

PS  ^  l/RS 


103 


^FA^AC  ~ 


FA^AC  ~  '^MAX,FA^AC  \j Xj PS  Xj RS 

X/AF  ^  X/PS  ^  X/RS 


0  4-1  yps 

YFA-^TG  '^MAX,FA^TG  y -\-X/PS^^^ 


10.  ~  ^1 


~ '^MAX,02^H^0  pc  7?C  <^2 

1  +  ^  +  ^ 

PS  RS 


yC^CR  ^MAX,PC^CR  p^  p^  ^ PC 


i  -A  r 

^CR^PC  ^^MAX,CR^PC  y p^  +XIPS^^^ 


Kg^fa  ^MAX,TG^FA^TG 


-k  _  2  '  r 

^LA^PY  ^MAX,LA^PY  yj^i^  X/RS 


j.  _  2  _ I _ p 

^CoApool^FC  ^MAX  ,CoApool^FC  Xj  RS  X/ FC 

l/RS  ^  l/FC 


16.  ^ AC^CoApool  ^MAX  ,AC^CoApool  Xj  PS  Xj  FC  ^ 

l/RS  ^  l/FC 


17.  <i>. 


-pisG 


ATP-yADP  ■  ^  R-ATP^ADP  1/ PS 
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